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SECTION  1 
INTRODUCTION 


This  report  presents  the  results  of  a  project  principally  aimed  at 
resolving  multivoice  data  into  separate  voice  signals.  A  secondary  aim 
was  to  resolve  noisy  multivoice  data  where  one  or  more  voices  may  be 
distorted.  Such  a  situation  might  result  from  adjacent  x^iio  channel 
interference  due  to  transmitter  or  receiver  nonlinearities . 

The  report  is  organized  mi inly  by  chronological  order.  Section  2 
describes  an  attempt  to  identify  to  which  of  two  speakers  various  sounds 
belong.  The  method  investigated  was  identification  by  adaptive  lattice 
filter  coefficients.  It  was  determined  that  this  is  not  practical. 

Section  3  describes  an  attempt  to  reduce  the  effect  of  a  primary 
speaker  masking  a  secondary  speaker.  The  method  used  was  'oeech  depend¬ 
ent  comb  filtering.  The  section  starts  with  a  description  of  complex 
correlation,  an  algorithm  used  extensively  throughout  the  successful 
phase's  of  the  project.  It  concludes  with  a  description  of  comb  filter¬ 
ing  and  its  partial  success. 

Section  4  returns  to  adaptive  lattice  filtering,  this  time  as  a 
possible  means  of  suppressing  a  primary  voice.  Listening  tests  showed 
that  it  did  not  separate  voices. 

Section  5  describes  the  sorting  of  multivoice  sounds  given  that 
certain  special  cases  are  satisfied.  Section  6  describes  the  combining 
of  the  work  of  Sections  3  and  5.  The  result  is  the  final  voice  processing 
developed  during  the  effort.  Section  6  also  describes  an  automatic  way 
to  estimate  speech  frequency  shift.  This  is  believed  to  be  a  new  result. 
Section  7  gives  conclusions,  and  recommendations  for  further  work. 
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SECTION  2 

STUDY  or  LPC  com' 'CIENTS 


The  first  stage  of  the  project  was  the  study  of  lattice  coefficients, 
also  known  as  linear  predictive  coding  (LPC)  coef f icients .  Other  terms 
for  them  are  reflection  coefficients  and  partial  correlation  (PARCOR) 
coefficients.  The  object  was  to  determine  if  the  coefficients  could  be 
used  to  assign  various  sounds  during  mixed  speech  to  their  respective 
speakers.  It  was  concluded  that  this  is  not  practical. 

2.1.  LATTICE  FILTERING 


The  starting  point  for  the  first  stage  of  the  project  was  a  paper 
by  Markel  et  al  [2]  entitled,  "Long-Term  Feature  Averaging  for  fpeaker 
Recognition."  In  this  paper,  the  authors  report  of  an  experiment  in 
identifying  speakers  by  pitch,  gain,  and  reflection  coefficients.  They 
found  the  not  surprising  result  that  the  longer  the  coef f icients  were 
averaged,  the  better  they  could  be  used  to  discriminate  between  speakers. 

The  authors  digitized  their  data  at  6.5  kHz,  and  applied  pre¬ 
emphasis,  before  performing  10  stage  lattice  filtering.  A  digitizing 
rate  of  6400  Hertz  was  selected  for  the  present  project.  Pre-emphasis 
was  applied  by  differencing  the  data  (1-z  J’)  before  lattice  filtering. 

The  number  of  samples  for  reflection  coefficient  computation  (the  frame 
size)  was  128  '20  milliseconds)  for  hoch  r he  first  stage  of  the  present 
project  and  the  work  of  Markel  et  al . 

The  method  used  in  the  present  project  for'  computing  reflection 
coefficients  was  Makhoul's  method  F  [1],  the  harmonic  mean  method, 
credited  to  Burg.  Makhoul  stated  that  ha  tended  to  prefer  the  use  of 
the  harmonic  mean  method  because  it  minimizes  a  reasonable  and  well- 
defined  error  criterion.  Figure  2'-l  illustrates  the  flow  of  data  through 
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Figure  2-1  Lattice 


a  lattice  filter  used  for  speech  analysis.  The  summer  at  the  left  is 
for  signal  pre-emphasir .  f^(n)  and  f^(n)  are  the  first  and  second 
forward  residuals  at  time  n.  b^(n)  and  b^Cn)  are  the  first  and  second 
backward  residuals  at  time  n.  k^  and  k^  are  the  first  and  second  reflec¬ 
tion  coefficients.  Coefficient  k^  at  stage  s  is  computed  as 

k  =  -2(Z  f  .  ( n  )b  (n-l)/(£  f2  (n)  +  I  b2  _  ( r— 1 ) ) 

5  S— 1  S-l  S-l  S-l 

n  n  n 

The  sum  over  n  was  taken  over  one  frame  interval,  namely  128  points  (20 

milliseconds),  for  this  stage  of  the  project.  The  reflection  coefficients 

are  calculated  in  order.  Once  reflection  coefficient  ko  has  been  computed, 

f  (n)  and  b  (n)  are  computed  for  all  n  in  the  frame  as 
s  s 

f£(n)  =  fs-l(n)  +  kSbS-l(n'1) 

and 

bs(n)  =  ksfs-l(n)  +  bs-l(n_1) 

Then  after  f  and  b  are  computed  k  may  be  computed.  As  Makhoul 

showed,  computing  as  just  given  minimizes  the  sum  of  the  squares  of 

f  and  b  . 
s  s 

Figure  2-2  shows  a  graphic  output  of  a  lattice  filtering  program 
applied  to  some  single  speaker  data.  The  lines  are  organized  in  groups 
of  three.  There  are  five  frames  per  line,  and  each  line  in  each  frame 
is  normalized.  The  first  line  of  each  group  shows  the  original  waveform 
data.  The  second  part  of  each  group  is  a  series  of  bar  graphs,  with 
each  bar  graph  showing  the  10  reflection  coefficients  for  or.e  frame. 

The  third  line  of  each  group  shows  the  forward  residual  after  10  stages 
of  lattice  filtering.  This  forward  residual  resembles  an  impulse  train 
during  voiced  speech. 
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2.2.  LATTICE  FILTERING  RESULTS 

Four'  male  speakers  were  recorded,  each  speaking  the  same  set  of  four 
sentences,  of  five  seconds  duration  each.  This  gave  a  total  of  16  data 
samples.  For  each  speaker  and  each  sentence,  a  two  second  average  was 
taken  for  each  of  the  ten  reflection  coefficients.  In  computing  each 
average...  a  weight  was  given  to  each  frame  proportional  to  the  signal 
power  in  the  frame.  This  was  to  reduce  the  weighting  of  reflection  coef¬ 
ficients  computed  during  silences  or  consonants. 

The  results  were  analyzed  with  the  On-Line  Pattern  Analysis  and 
Recognition  System  (OLPARS).  It  was  found  that  coefficients  3  and  6  were 
best  for  discriminating  between  speakers.  This  is  similar  to  the  results 
of  Markel  et  al,  who  found  that  coefficients  2  and  6  gave  the  best 
discrimination.  The  difference  in  result  may  be  due  to  the  small  data 
base  of  the  present  effort.  Figure  2-3  shows  an  OLPARS  plot  in  the 
plane  of  reflection  coefficients  3  and  6.  Letters  A  through  D  represent 
the  four  speakers.  The  four  data  points  for  each  speaker  are  connected 
by  a  polygon  drawn  by  hand. 

The  OLPARS  plot  shows  that  in  this  case  speakers  A,  B,  and  D  can  be 
distinguished,  while  speaker  C  data  overlaps  the  data  for  both  speakers 
A  and  D.  This  result  is  not  encouraging,  because  two  seconds  of  data  is 
still  much  more  than  can  be  obtained  from  a  single  syllable.  Markel  et 
al  determined  that  reflection  coefficient  scattering  is  inversely  propor¬ 
tional  to  the  cube  root  of  the  number  of  voiced  samples  averaged. 
Therefore,  decreasing  the  averaging  interval  by  a  factor  of  27  (from  2 
seconds  to  74  milliseconds)  would  increase  the  scatter  for  each  speaker 
by  a  factor  of  three.  This  would  cause  the  data  for  the  four  speakers  to 
overlap  substantially. 
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An  additional  problem  is  that  reflection  coefficients  are  sensitive 
to  noise  and  to  competing  speaker  signals.  It  was  concluded  that  sorting 
syllables  based  on  reflection  coefficients  would  not  be  reliable.  Addi¬ 
tional  work  with  LPC  lattice  filtering  is  discussed  in  Section  u . 
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SECTION  3 

SPEECH  DEPENDENT  COMB  FILTERING 


3.1.  COMPLEX  CORRELATION 

By  complex  correlation  is  meant  an  operation  similar,  but  noi. 
identical,  tc  computing  an  autocorrelation.  An  autocorrelation  may  be 
computed  by  taking  an  FFT,  followed  by  taking  the  magnitude  squared, 
followed  by  taking  an  inverse  FFT.  As  the  term  will  be  used  here, 
complex  correlation  starts  with  taking  a  tapered  window  of  the  data, 
namely  the  Hanning  (raised  cosine)  window.  This  is  followed  by  taking 
an  FFT,  and  then  taking  the  magnitude  squared.  The  result  is  the  power 
spectrum.  This  corresponds  so  far  to  the  operations  in  calculating  an 
autocorrelation.  Next,  however,  for  taking  a  complex  correlation,  the 
negative  frequency  portion  of  the  transform  is  zeroed.  For  the  discrete 
Fourier  transform  of  N  points,  this  is  equivalent  to  zeroing  the  values 
at  N/2,  N/2+1,  ...,  N-l .  No  information  is  lost  by  this  operation, 
because  the  power  spectrum  of  a  real  signal  is  symmetric  about  zero 
frequency.  The  discrete  power  spectrum  of  a  pure-real  sampled  time 
signal  is  also  symmetric  about  half  the  sampling  frequency.  This  means 
that  the  discrete  power  spectrum  values  at  N/2,  N/2+1,  ...,  N-l  are  the 
mirror  images  of  the  values  at  0,  1,  ...,  N/2-1.  Thus  when  half  the 
power  spectrum  is  zeroed  only  redundant  information  is  lost. 

Next  in  computing  the  complex  correlation,  the  square  roots  of  the 
remaining  points  are  taken.  This  operation  is  done  to  reduce  the  ratio 
of  peaks  in  the  data.  For  example,  a  power  spectrum  with  one  peak  25 
times  as  high  as  another  will  have  this  ratio  reduced  to  5  after  the 
square  root  is  taken.  Since  the  complex  correlation  will  be  used  to 
estimate  pitch,  it  is  important  that  it  be  influenced  by  several  pitch 
harmonics,  and  not  just  by  a  single  dominant  power  spectrum  peak. 
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Finally,  the  inverse  FFT  is  taken.  The  result  is  a  complex  func¬ 
tion,  because  it  is  the  inverse  transform  of  a  function  with  only  posi¬ 
tive  frequency  components.  Hence  the  term  complex  correlation.  One  way 
of  viewing  the  complex  correlation  is  that  each  point  in  it  represents  a 
coefficient  in  the  Fourier  series  which  can  reproduce  the  processed 
power  spectrum.  For  S(0),  S(l),  ...,  S(N-i)  denoting  the  processed 
power  spectrum,  we  have  that  the  correlation  Ct.0),  C(l),  . ..,  C(N-l)  is 
given  by 

1  N-X 

C(n)  =  —  Z  S(k)  exp(  j2trkn/N) 
k=0 

for  n  =  0 ,  1 ,  . . . ,  N-l . 

See  [3]  page  89.  Then  S(k)  is  represented  by  the  Fourier  series 
with  coefficients  C(n),  that  is  by 

N-l 

S(k)  -  T.  C(n)  exp(- j 27Tkn/N) 
n  =  0 

However,  S(k)  is  pure  real,  so  we  may  rewrite  this  as 
N-l 

S(k)  =  z  Re(C(n))  cos(2iTkn/N) 
n=0 

N-l 

+  £  Im(C(n))  sin(2irkn/N), 

n=0 

3.2.  COMPLEX  CORRELATION  APPLIED  TC  SPEECH  DATA 


Speech  data  during  voiced  speech  is  characterized  by  the  presence 
of  multiple,  harmonics  of  the  pitch.  During  unvoiced  speech,,  the  data 
is  characterized  as  filtered  broadband  noise.  However,  the  unvoiced  por¬ 
tions  are  relatively  less  important  for  two  reasons.  First,  unvoiced 
speech  has  considerably  less  power  than  does  voiced  speech.  Thus  in 
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noisy  data,  it  is  to  a  considerable  extent  masked  by  the  noise.  Second, 
listening  tests  have  shown  that  the  human  ear  is  relatively  insensitive 
to  distortion  in  the  reproduction  of  unvoiced  sounds.  For  these  reasons, 
the  speech  processing  under  the  present  effort  was  limited  to  treating 
all  data  as  if  it  consisted  of  voiced  speech.  Listening  tests  of  the 
results  confirm  that  the  resulting  suboptimal  processing  of  unvoiced 
sounds  is  hardly  noticeable.  Indeed,  it  seems  likely  that  attempting  to 
make  voiced/unvoiced  decisions  on  noisy  multitalker  data  would  introduce 
more  distortion  (due  to  decision  errors)  than  it  would  eliminate. 

Figure  3-1  illustrates  the  square  root  of  power  spectrum  versus 
time  for  a  single  talker.  This  is  for  microphon  :  speech  (as  opposed  to 
SSB  radio  reception)  with  only  ambient  room  noise  in  the  background. 

Each  line  represents  a  frequency  range  of  0  to  3200  Hertz,  and  time 
advances  by  40  milliseconds  from  one  line  to  the  next.  Many  lines  show 
the  comb-like  regular  succession  of  peaks  characteristic  of  voiced 
speech. 

Figure  3-2  illustrates  the  magnitude  of  a  complex  correlation 
versus  time  for  a  single  talker.  This  is  also  for  microphone  speech. 

This  plot  was  produced  by  the  program  CCORPL,  listed  in  Appendix  B. 

Each  horizontal  line  represents  a  time  range  of  0  to  20  milliseconds. 

Time  advances  by  20  milliseconds  from  one  line  to  the  next.  The  curve 
to  the  left  of  the  horizontal  lines  shows  speech  amplitude  versus  time. 
Correlation  peaks  are  visible  on  most  of  the  horizontal  lines,  showing 
pitch  periods  of  8  to  10  milliseconds.  The  speaker's  pitch  is  therefore 
in  the  range  of  100  to  125  Hertz.  Several  portions  of  time  where  the 
speaker's  pitch  is  not  evident  also  have  dips  in  speech  amplitude. 
Absence  of  correlatio1-  peaks  showing  speaker  pitch  indicates  silence,  or 
non-voiced  speech.  The  complex  correlation  provides  an  automated  way  of 
determining  speaker  pitch.  It  may  also  be  used  to  determine  the  best 
fit  of  a  comb  filter  to  the  data.  This  will  be  further  discussed  in  the 
next  subsection. 
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3.3. 


PRINCIPLES  OP  SPEECH  DEPENDENT  COrtB  FILTERING 


By  speech  dependent  comb  filtering  is  meant  applying  a  filter  to 
the  data  with  the  filter  itself  dependent  on  the  data.  In  particular, 
the  object  is  to  enhance  or  suppress  the  voiced  speech  oi  one  speaker. 

As  Figure  3-1  illustrates,  voiced  speech  has  a  comb-like  structure  in 
the  frequency  domain.  To  enhance  such  data,  a  filter  may  be  used  which 
passes  the  speech  peaks  and  attenuates  the  valleys.  To  suppress  the 
speech,  a  filter  may  be  used  which  attenuates  the  peaks  while  parsing 
the  valleys.  The  purpose  of  suppressing  the  speech  is  to  let  the  speech 
of  a  second  speaker  come  through  with  comparatively  little  attenuation 
versus  the  primary  speaker. 

From  Section  3.1,  we  have  that  for  S(k)  the  one-sided  square  root 
of  the  power  spectrum  and  C(n)  f  . omplex  correlation 

L  N-l 

-  N"  ^  S(k)  exp( jGirkn/N) 
k=0 

We  wish  to  find  the  pitch  period  T  and  frequency  displacement  d 
such  that  the  comb  function  F(k)  best  fits  S(k)  where 

F(k)  -  0.5  +  0.5  cos( 2irTk  +  d) 

That  is,  we  wish  to  find  the  T  and  d  which  maximize  the  dot  product  of 
F(k)  and  S(k): 

N-l 

F  •  S  =  £  0.5[1  +  Re(exp(  j2TrTk)exp(jd)  )]S(k) 

k=0 


-  0. 5NC( 0 )  +  0. 5NRe(C(T)  exp(jd)) 
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The  maximizing  T  may  not  be  an  integer.  We  will  deal  with  this  problem 
shortly.  Given  'I  and  a  complex  C(T),  then  the  maximizing  d  is  any  value 
which  makes  C(T)  exp(jd)  pure  real.  It  follows  that 

max  F*S  =  0.5N(C(0)  +  |c(T) | ). 
d 

Therefore  the  T  which  gives  the  best  comb  fit  to  the  data  is  the  T  which 
gives  the  maximum  value  of  the  magnitude  of  the  complex  correlation.  Of 
course,  we  must  rule  out  very  small  values  of  T,  because  the  complex  cor¬ 
relation  C(n)  has  its  overall  maximum  value  at  n=0.  We  are  interested 
only  in  the  maximum  value  over  the  range  of  possible  pitch  periods. 

The  d  which  gives  the  best  fit  is  such  that  exp(jd)  =  C(T)!’!/  j  C(T)  |  , 
where  ■'  denotes  the  comp1  ex  conjugate.  That  is, 

cos  d  -  Re  C(T )/ jC(T) j  and 
sin  d  =  -im  C(T )/ | C(T)  | 

It  follows  that 

d  =  arctan(-Im  C(T)/R3  C(T)). 

Care  must  be  taken  in  evaluating  the  arctan  to  put  d  into  the  proper 
quadrant.  In  Fortran,  this  may  be  done  by  using  the  ATAN2  function. 

Now  we  consider  the  choice  of  a  non-integer  T.  The  first  method 
used  to  select  the  best  T  was  to  first  find  the  maximizing  integer  T^, 
and  then  take  the  centroid  (center  of  gravity)  of  the  three  points  C(T^ 

-  1),  C(T^)  and  C(I\  +  1).  This  gave  the  centroid  value  T£  as 

(T.  -  l)C(T .  -  1)  +  T.C(T.)  +  (T.  +  1  )C(T.  +  .1) 

T  _  __i _ i _ ii  i _ i 

c  "  C(T.  -  1)  +  C(T . )  +  C(T .  +  1)  * 
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Then 


C(T .  +  1)  -  C(T .  -  1) 

T  _  T  _  i  i  _ 

'c  i  "  C(!T.  -  l)  +  C(T.  )  +  c<!t .  +  l)  ' 
i  ii 

The  method  used  to  find  the  best  d  given  a  non-integer  T  was 
quadratic  interpolation.  Three  points  determine  a  quadratic  curve,  and 
the  three  points  used  were  the  best  integer  T  (denoted  T^),  along  with 
Tj.  -  1  and  T.  +  1.  The  real  and  imaginary  parts  of  C(T)  were  interpo¬ 
lated  separately. 

For  G(k)  a  function  of  integer  valued  k,  we  may  fit  to  it  the 

2 

quadratic  function  F(x)  =  G(0)  (-  C^x  +  C2X^  w^2re 
C,  =  0. 5(G( 1)  -  G (-1)) 

C2  =  0. 5(G(-l)  +  G(l) )  -  G(0 ) 

The  reader  may  verify  that  F(x)  =  G(x)  for  x  =  -1,  0,  1 .  Then  for 
T^  the  best  integer  T,  and  Tfa  the  best  (possibly  non-integer)  T,  the 
real  part  of  C(Tfa)  was  taken  to  be 

Re  C(T.)  =  C(T.)  +  (T,  -  T.)  ( R  +  (?.  -  T.)R„) 
b  l  bil  Dir 

where 

R.  =  0.5  Re(C(T.  +  1)  -  C(T.  -  1)) 

1  ii 

R„  =  Re[0.5(C(T.  -  1)  +  C(T.  +  1))  -  C(T.)]  . 

^  ill 

The  imaginary  part  of  C(T^)  was  computed  similarly.  These  methods  of 
interpolation  were  used  through  the  time  that  the  first  few  automatic 
mistuning  estimates  were  made  (see  Section  6.1).  It  was  found  then  that 
the  mistuning  estimates  were  consistently  about  -5  Hertz  for  microphone 
speech,  which  had  never  been  mistuned.  A  revision  was  made  to  the 


3-8 


r 


method  of  computing  the  optimal  pitch  period  T.  The  revised  idea  was  to 

choose  the  T  which  maximized  the  quadratic  interpolation  of  |C(n)| 

fitted  to  the  three  points  |C(Ti  -  D  J C( )  (  ,  and  |C(Ti  +  1^2‘ 

This  T  may  be  determined  as  follows.  For  F(x)  =  G(0)  +  C^x  +  C?x  we 
o 

have 


F'(x)  =  C1  +  2C2x 


F" (x)  =  2C, 


so  F(x)  has  its  maximum  at  x  =  -C1/2C2>  provided  that  C2  is  negative. 
For  the  quadratic  interpolation  of  |C(n)j  we  have 


x  =  T  -  T. 

o  i 


Cx  =  0.5(|C(T.  +  1)|2  -  |C(T.  -  1)|2) 

C  =  0.5(|C(Ti  -  l)i2  +  (C(T i  +  1)(2)  -  (C(Ti) j2 


Therefore  we  choose 


T  =  T.  ~  C./2C_ 

o  1  1  2 


provided  that  |c(T.)|2  >  0.5(|C(T.  -  l)|2  +  |C(T.  +  1)|2)-  This  last 

condition  is  true  any  time  it  is  true  that  |C(T.)r  is  greater  than  both 

|C(T.  -  1)|2  and  lc^Ti  +  L^|2-  UsinS  the  quadratic  interpolation  value 

for  T  was  found  to  correct  the  mistuning  estimates, 
o 
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VARYING  COMB  TOOTH  THICKNESS 


the  resulting  signal  still  had  its  spectral  peaks  where  the  original 
signal  had  them.  This  was  believed  to  t>e  due  to  frequency  quantizing 
preventing  the  voice  suppression  filter  (referred  to  as  a  rake  filter) 
from  ever  being  exactly  zero.  Figure  3-3  illustrates  the  phenomenon. 

Possibly  widening  the  valleys  of  the  rake  filter  would  suppress  the 
voiced  speech  even  more  effectively,  and  as  a  result,  let  some  secondary 
speech  come  through  better.  Similarly,  if  a  secondary  speech  peak 
coincided  with  a  rake  filter  peak,  then  possibly  widening  the  rake 
filter  peak  (its  "tooth")  would  allow  the  secondary  speech  signal  to 
come  through  better. 

The  method  used  to  make  peaks  ana  valleys  wider  or  narrower  was  to 
warp  the  frequency  axis  when  applying  ■•'he  comb.  That  is,  instead  of  the 
comb 


F(k)  =  0.5  +  0.5  cos(2nTk  +  d) 
the  frequency  warped  comb  F(z(k))  was  used: 

F( z(k ) )  =  0.5  *  0.5  cost 2mTz(k)  +  d). 

Warpad  frequency  z(k)  was  determined  as  follows:  For  each  frequency 
point  k,  the  points  x^  and  x^  are  such  that  x^  is  the  location  of  the 
comb  peek  or  valley  nearest  to  but  preceding  k.  x^  is  the  peak  or 
valley  nearest  to  but  following  k.  Then  if  xQ  is  a  peak  location,  x^  is 
a  valley  location,  and  vice  versa.  Then  for  S(n)  the  square  root  of  the 
power  spectrum,  z(k),  is  computed  to  be 

z(k)  =  k  +  K(S(xQ)  -  S(x  ))(k  -  xQ)(k  -  x^ 
where  K  is  such  that 
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1/K  =  (S(x0)  +  S(x  iHx^  -  >'0 ) 


For  non- integer  Xq  and  x^  linear  interpolation  is  used  to  compute  S(Xq) 
and  S(x^)  respectively.  Then  z(k)  =  k  for  k  =  x^  or  x^. 

A1S°  xQ  +  Xj  xQ  +  x^  (S(x^)  -  S(xQ))(x1  -  xQ) 

z(  2  }  =  2  +  4^S(xQ)  +  S(xx)) 

xo  +  xi  x0  +  X1 

So  z( - - — — )  >  — -  if  and  only  if  Stx^)  >  S(xQ). 


Thus  if  a  power  spectrum  peak  matches  either  a  comb  peak  or  valley,  then 
that  comb  peak  or  valley  is  widened  to  include  more  of  the  power  spectrum 
peak.  The  coefficient  K  in  the  warping  formula  is  the  maximum  value 
consistent  with  not  warping  points  in  the  interval  xQ  to  x  so  much  that 
they  are  moved  outside  the  interval.  This  is  equivalent  to  the  condition 
that  z' (xQ)  >_  0,  and  z’(x^)  _>  0.  We  have 


Now 


z'(k)  =  1  +  K(S(xq)  -  S  ( x ,  ) )  ( k  -  xQ  +  k  -  x^ 

S(xQ)  -  S(x1) 

_1  -  S(xQ)  +  S(xx)  -  1 


so 


and 


z  (x  )  >  1  +  ( - )  (x  -  x,  )  =  0 

0  —  x.  -  x„  0  1 

x  0 


z'(x  )  >  1  -  ( - )  (x  -  x  )  =  0. 

1  —  x1  -  x0  1  0 


3.5.  COMB  VOICE  PROCESSOR  AND  RESULTS 


Figure  3-4  gives  a  block  diagram  of  the  data  flow  within  the  comb 
voice  processor.  The  input  data  is  digitized  at  6400  Hertz.  It  is 
divided  into  512  point  (80  mio-lisecond )  intervals,  with  each  interval 


advanced  256  points  from  the  previous  one.  That  is,  the  intervals  have 
50  percent  overlap.  Each  window  is  given  the  Hanning  (raised  cosine) 
weighting,  and  a  fast  Fourier  transform  (EFT)  is  taken.  This  transform 
is  saved  for  further  use.  The  square  root  of  the  power  spectrum  is 
computed  from  the  EFT  with  negative  frequency  portions  zeroed.  Next  an 
inverse  FFT  is  taken  of  the  root  power  spectrum.  This  is  the  complex 
correlation  of  the  signal,  as  described  in  Section  3.1.  Next  the  magni¬ 
tude  peak  is  found  for  the  complex  correlation ,  and  the  phase  at  the 
peak  is  computed.  The  search  is  performed  over  points  25  to  80  of  the 
complex  correlation,  corresponding  to  speaker  pitches  of  81  to  267 
Hertz.  The  location  of  the  correlation  peak  and  the  phase  at  the  peak 
allow  determination  of  the  best  comb  filter  fit  to  the  data,  as  described 
in  Section  3.3. 

At  this  point,  frequency  warping  may  be  applied.  A  comparison  is 
made  of  the  linear  magnitude  power  spectrum  versus  the  locations  of  the 
comb  filter  peaks  and  valleys.  The  frequency  axis  is  warped  accordingly, 
as  described  in  Section  3.4.  The  same  frequency  warping  is  used  both 
for  speech  enhancement  and  suppression.  This  is  because  the  enhancement 
comb  peaks  line  up  with  the  suppression  comb  valleys,  and  the  frequency 
warping  does  not  distinguish  between  peaks  and  valleys. 

The  FFT  of  the  original  signal  is  multiplied  by  the  best  fit  raised 
sinusoid  (possibly  with  frequency  warping)  and  the  inverse  FFT  is  taken 
to  give  the  primary  output.  Multiplication  in  the  frequency  domain  is 
equivalent  to  convolution  in  the  time  domain,  so  the  effect  is  a  linear 
filtering  of  the  signal.  A  secondary  output  is  produced  by  multiplying 
the  signal  FFT  by  a  raised  sinusoid  (possibly  with  frequency  warping) 
with  the  same  spacing  between  the  teeth  as  for  the  primary  output. 
However,  for  the  secondary  output,  the  comb  is  positioned  so  that  overall 
its  teeth  line  up  against  valleys  in  the  signal. 
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The  comb  voice  processor  by  itself  was  a  partial  success,  as 
Figures  3-5  through  3-9  illustrate.  Each  of  these  figures  shows  the 
evolution  of  a  0  to  3200  Hertz  spectrum  versus  time.  Each  line  is 
derived  from  a  40  millisecond  (256  point)  Hanning,  or  raised  cosine, 
time  window.  There  is  a  10  millisecond  advance  between  lines.  Figures 
3-5  and  3-6  show  spectra  for  speakers  A  and  B,  respectively.  Figure  3-7 
shows  spectra  for  the  speech  of  speaker  A  mixed  at  6  dB  above  (at  four 
times  the  power  of)  the  speech  of  speaker  B.  The  B  speech  shows  up 
mainly  as  a  perturbation  to  the  A  speech.  Figure  3-8  shows  spectra  for 
the  mixed  speech  subjected  to  the  suppression  algorithm  without  frequency 
warping.  Figure  3-9  is  for  suppression  with  frequency  warping.  The 
result  of  suppression  without  frequency  warping  shows  considerably  more 
evidence  of  speaker  B  than  does  the  original  mixed  speech.  Comb  filtering 
is  there  ore  at  least  a  partial  success.  The  result  of  suppression  with 
frequency  warping  is  inferior  to  the  result  of  suppression  without 
frequency  warping.  For  this  reason,  frequency  warping  wa^  not  used  in 
the  later  stages  of  the  project.  Incidentally,  the  voice  processing 
whose  results  are  shown  here  used  the  quadratic  interpolation  method  of 
calculating  speaker  pitch. 

Also  tried  was  saturating  (clipping)  the  raised  sinusoid  frequency 
response  so  that  it  is  zero  at  25%  of  the  frequencies,  one  at  25%  of  the 
frequencies,  and  transiting  between  zero  and  one  for  50%  of  the  frequen¬ 
cies.  There  seemed  to  be  a  slight  but  not  significant  degradation 
versus  the  raised  sinusoid  comb.  The  saturated  comb  was  not  used  in  the 
later  stages  of  the  project. 

The  initial  hope  for  the  comb  voice  processor  was  that  it  would  be 
enough  by  itself  to  make  intelligible  a  secondary  voice,  which  prior  to 
processing  had  been  masked  by  a  primary  voice.  This  turned  out  not  to 
be  the  case.  The  probable  cause  of  the  failure  is  that  there  are  signi¬ 
ficant  amounts  of  time,  perhaps  20  to  30  percent,  when  the  primary 
speaker  is  not  producing  sound,  but  the  secondary  speaker  is.  During 
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Suppression  with  Frequency  Warping 


these  times  a  pure  suppression  filter  will  lock  on  to  the  secondary 
voice  and  suppress  it.  What  is  needed  is  a  way  for  the  signal  pro¬ 
cessing  algorithm  to  distinguish  intervals  when  one  voice  is  dominant, 
versus  intervals  when  the  other  voice  is  dominant.  Two  methods  for 
doing  this  for  special  cases  are  discussed  in  Sections  5  and  6. 
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SECTION  4 


LPC  ANALYSIS  AND  RECONSTRUCT i ON 


Section  2  discussed  LPC  lattice  filtering  with  the  object  of  assign¬ 
ing  syllables  in  mixed  speech  data  to  their  respective  speakers  .  Here  we 
discuss  the  use  of  lattice  filtering  as  a  stage  of  speech  processing. 

In  this  portion  of  the  project,  lattice  filtering  was  followed  by 
thresholding  the  forward  residual,  followed  by  waveform  reconstruction . 
The  hope  was  that  the  thresholding  would  pass  a  dominant  voice  while 
rejecting  a  secondary  voice.  The  reconstructed  signal  might  then  be 
subtracted  from  the  original  data,  leaving  trie  secondary  voi<ce.  It 
turned  out  that  the  procedure  did  not  separate  the  voices. 

4.1.  SPEECH  RECONSTRUCTION  3Y  LATTICE  FILTERING 


Figure  2-1  illustrates  the  initial  stages  of  a  lattice  analysis 
filter.  As  described  in  Section  2.1,  the  forward  residual  f  is  related 
to  the  backward  residual  b  by 


and 


f  (n)  =  f  (n)  +  k  b  ,  ( n  - 1 ) 
S  s-l  S  S-l 


b  ( n )  -  k  f  In)  +  b  , ( n- 1 ) . 
s  s  s-l  s-l 


We  may  solve  these  equations  for  (n)  and  b_(:i>,  giving 


f  ,(n)  -  f  '  n.  -  k  b  .  (n-1) 
s-I  s  c  s-i 


and 


b  (n)  -  k  f  (n)  +  ( ' -k  b  .(n-1), 
s  s  s  s  s-l 


Figure  4-1  illustrates  how  these  equations  may  be  implemented.  The 
processing  of  Figure  4-1  is  the  inverse  of  the  processing  of  Figure  2-1. 
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The  summer  circuit  on  the  left  in  4-1  is  the  reciprocal  of  the  differencing 
circuit  on  the  left  in  Figure  2-1.  The  speech  synthesis  circuit  of 
Figure  4-1  requires  only  the  forward  residual  from  the  analysis  filter 
of  Figure  2-1. 

In  the  analysis  filter,  the  s-th  residuals  fg(n)  and  b^tn)  are  com¬ 
puted  for  fixed  s  and  for  all  times  n  before  the  computation  of  fg  (n) 

and  b  ,(n).  By  contrast,  in  the  synthesis  filter,  f  (n)  and  b  (n)  are 
s+1  s  s 

computed  for  fired  time  n  and  all  stages  s  before  the  computation  of 

f  (n+1)  and  b  (n+1).  The  synthesis  filter  memories  are  assumed  to  be 
s  s 

filled  with  zeroes  initially. 

4.2.  LPC  VOICE  PROCESSOR  £ND  RESULTS 


The  voice  analysis  and  reconstruction  described  in  Section  4.1  was 
programmed.  Twenty  stages  of  lattice  filtering  were  used.  The  theory 
was  that  ten  stages  are  suitable  for  a  one  voice  signal,  and  a  two  voice 
signal  may  be  Twice  as  complex  as  a  one  voice  signal.  Data  segments  of 
149  points  (23.1  milliseconds)  were  used,  with  20  point  (3.1  millisecond) 
overlaps  between  segments,  to  allow  for  lattice  filter  startup  effects. 
The  analysis  filtering  of  each  segment  resulted  in  a  20  point  startup 
ara  a  128  point  forward  residual. 

The  root  mean  square  (rms)  value  of  the  forward  residual  was  calcu¬ 
lated.  Then  the  magnitude  of  each  point  in  the  residual  was  compared  to 
a  threshold  consisting  of  the  rms  value  times  a  constant.  Primary  and 
secondary  I'esiduals  were  then  computed.  The  primary  residual  contained 
all  residual  values  whose  magnitudes  exceeded  the  threshold,  and  was 
filled  with  zeroes  elsewhere.  The  secondary  residual  contained  all 
residual  values  whose  magnitudes  fell  below  thr  threshold,  and  was 
filled  with  zeroes  elsewhere. 


4-3 


Figure  4-2  shows  a  result  from  the  processing  of  a  two-voice  signal 
with  this  algorithm.  In  this  case,  the  threshold  for  the  magnitude  of 
the  residual  was  set  at  one  times  the  rms  value.  The  lines  of  the  page 
are  organized  into  groups  of  four.  There  are  five  frames  per  line,  and 
each  line  in  each  frame  is  normalized.  The  first  line  of  each  group  shows 
the  original  waveform  data.  The  second  line  of  each  group  shows  the 
forward  residual  after  20  stages  of  lattice  filtering.  The  third  line 
shows  the  primary  residual  after  thresholding.  The  fourth  and  last  line 
of  each  group  shows  the  primary  output,  that  is,  the  output  synthesized 
from  the  primary  residual. 

A  program  was  written  which  produced  primary  and  secondary  output 
files.  This  program  was  named  VCRFSP  and  is  listed  in  Appendix  B.  The 
program  was  applied  to  mixed  speaker  data.  Listening  tests  on  the  results 
showed  that  the  primary  output  generally  contained  both  voices.  The 
secondary  output  contained  both  voices,  but  sounded  as  if  both  talkers  were 
whispering.  There  is  a  logical  explanation  for  this  phenomenon. 

A  person  speaking  voiced  speech  produces  an  impulse  train  at  his 
vocal  cords,  and  this  impulse  train  acoustically  excites  his  vocal  tract. 

By  contrast,  a  whisperer  does  not  produce  impulses  at  his  vocal  cords. 
Instead  he  excites  his  vocal  tract  with  the  white  noise  generated  by  the 
rush  of  air  through  his  throat.  Of  course,  a  voicing  person  also  excites 
his  vocal  tract  with  white  noise,  but  the  effect  is  masked  by  the  higher 
amplitude  excitation  of  the  impulse  train.  Lattice  filtering  in  effect 
undoes  the  acoustic  filtering  of  the  vocal  tract,  resulting  in  a  forward 
residual  representing  the  excitation  of  the  vocal  tract.  Then  removing 
the  impulses  from  this  excitation  will  remove  the  voiced  portion  of  the 
signal.  Reconstruction  will  then  result  in  making  audible  the  whispering 
portion  of  a  voicing  speaker's  sound.  Unfortunately,  this  does  not  help 
separate  voices . 
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A  possible  further  course  of  research  along  LPC  lines  remains, 
is  to  try  to  sort  out  mixed  impulse  trains  in  the  forward  residual, 
assigning  one  series  of  impulses  to  one  speaker,  and  another  series  to 
another  speake  .  This  line  of  research  was  not  pursued  for  three 
reasons.  First,  sorting  impulses  would  be  difficult  because  the  impulses 
are  not  always  clearly  distinct  from  the  low  level  noise.  Second,  even 
if  the  impulses  could  be  isolated,  and  impulse  trains  separated,  there 
would  be  the  problem  of  identifying  which  speaker  produced  which  impulse 
train.  Third,  the  reflection  coefficients  are  affected  by  both  vocal 
tracts.  Even  if  the  impulse  trains  could  be  properly  assigned,  there 
would  be  the  problem  of  resolving  the  lattice  filter  coefficients  in 
some  way,  so  as  to  reconstruct  each  speaker’s  speech  separately,  without 
influence  from  the  shape  of  the  competing  speaker's  vocal  tract.  Because 
of  these  difficulties,  further  LPC  research  was  not  pursued  during  the 
project. 
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SECTION  5 
SYLLABLE  SORTING 


The  idea  of  syllable  sorting  is  to  assign  each  section  of  mixed 
voice  data  to  one  or  the  other  of  two  output  channels,  depending  on 
which  speaker's  voice  is  dominant  in  the  data  section.  No  attempt  is 
made  to  split  a  sound  into  parts.  Listening  tests  showed  that  the 
method  can  separate  mixed  voices  fairly  well  if  a  good  sorting  criterion 
is  available.  Such  criteria  were  developed  for  two  special  cases.  The 
resulting  speech  processing  is  discussed  in  sections  5.1  and  5,2. 

5.1.  SORT  BY  PITCH  ALGORITHM 

One  obvious  special  case  which  allows  assigning  sounds  to  speakers 
is  one  where  the  speakers  have  distinct  pitch  ranges.  This  might  typi¬ 
cally  be  found  in  mixed  speech  from  a  male  and  a  female  speaker.  For 
each  section  of  data,  the  predominant  pitch  may  be  determined  by  complex 
correlation  with  a  search  for  the  peak,  as  described  in  sections  3.1  and 

3.2.  One  advantage  of  using  the  complex  correlation  is  that  it  is 
robust  against  frequency  shift  of  the  signal,  such  as  could  result  from 
single  sideband  (SSB)  m.istuning. 

Sorting  by  pitch  was  implemented  using  overlapping  Hanning  (raised 
cosine)  time  windows  applied  to  the  data.  The  sampling  rate  was  6400 
Hertz  and  the  window  overlap  was  50  percent.  Window  lengths  of  256  and 
512  points  (40  and  80  milliseconds)  were  tried.  Each  window  was  assigned 
to  one  of  two  output  channels,  depending  on  the  window's  dominant  pitch. 
The  output  channel  which  was  rejected  received  a  section  of  white  noise, 
at  about  6  dB  below  the  mixed  speech  "ignal.  This  white  noise  was 
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inserted  to  prevent  gaps  in  the  output  channels.  It  is  difficult  for  a 
human  to  interpret  speech  interrupted  by  silences. 

The  algorithm  separated  speech  well,  provided  that  the  speakers  had 
distinct  pitch  ranges.  Use  of  the  80  millisecond  windows  was  found  to  be 
superior,  because  occasionally  the  algorithm  would  switch  back  and  forth 
rapidly  between  the  output  channels.  The  algorithm  with  the  40  millisecond 
windows  would  then  produce  a  rapid  (50  Hertz)  fluttering  effect,  which  was 
annoying  and  distracting. 

To  determine  how  well  the  windows  were  being  sorted,  an  ideal  syllable 
sorter  was  programmed.  This  ideal  sorter  was  given  two  separate  voice 
files,  which  it  then  merged  at  a  specified  power  ratio,  ’’’hat  is,  th 
program  was  told  how  many  dB  the  first  voice  should  be  above  the  second  in 
the  mixed  file.  The  program  then  computed  the  mixed  voice  signal  and  com¬ 
pared  it  versus  the  two  separate  voice  files.  For  each  window  of  data, 
the  program  computed  the  gain  which  would  minimize  the  mean  square  error 
(mse)  between  the  separate  voice  file  and  the  mixed  voice  file.  The 
program  then  computed  mse's  and  assigned  the  data  accordingly. 

The  optimal  gain  to  minimize  mean  square  error  may  be  computed  as 
follows.  For  two  real  functions  f  and  g  define  the  inner  product  of  f 
and  g  as 

N 

<f,g>  =  E  f(n)  g(n)  0.5(1  -  cos(27rn/N)) 
n=l 

Define  the  norm  |  |f|  |  =  Af  ,f  > 

It  may  be  verified  that  these  are  a  true  inner  product  and  norm.  Then 
given  a  solo  voice  signal  f  and  a  mixed  voice  signal  g,  we  wish  to  find 
the  gain  constant  a  which  minimizes  the  mse,  given  by 
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mse  =  |  |af  -  p  |  |A  -  <af~g,  af~g>  - 

Here  we  make  use  of  the  property  that  <f,g> 
Then 


d(mse) 

da 


2a|  |f j  |2  -  2 <f , g > 


and 


d  (mse) 
da2 


=  2  f 


!  U 


=  <g ,  f  >. 


2a<f,g>  t 


The  second  derivative  is  non-negative,  so  the  mse  has  its  minimum  where 

its  first  derivative  equals  zero.  Let  a  be  the  a  which  minimizes  the 

o 

mse .  Then 


a 

o 


<f.g>/| lf I  I' 


and  the  minimum  mse  is 


min(mse)  =  |jg||2  -  <f  ,g>2/' |  |f  |  j  2 . 
a 


The  power  of  the  modified  solo  voice  signal  a  f  is  a 

o  o 

best-.fit  signal  to  noise  ratio  (SNR)  is 


so  the 


SNR 


|2  -  <f,g>2/||f| 


So 


SNR  =  - 

ilfll'ligli' 

provided  ||f||2||g||2  >  0. 


The  ideal  syllable  sorter  was  given  a  tolerance  in  dfl .  It  sorted 
each  mixed  voice  window  into  the  output  channel  representing  the  speaker 
with  the  higher  SNR  for  the  window.  It  assigned  the  mixed  voice  signal 
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also  to  the  alternate  output  channel  if  both  SNR's  were  within  the  toler¬ 
ance  of  each  ether.  If  the  SNR's  differed  by  more  than  the  tolerance,  then 
the  alternate  output  channel  was  given  white  noise. 

Tor  voices  mixed  at  0  dB  relative  levels,  it  was  found  that  ideal 
sorting  with  10  dB  tolerance  allowed  nearly  complete  reconstruction  of 
each  voice,  with  only  a  small  amount  of  the  alternate  voice  coming  through 
on  each  channel.  Thus  syllable  sorting  alone  is  able  to  recover  all  of  a 
person's  speech  that  is  no  more  than  10  dB  below  a  competing  talker's 
speech.  This  implies  that  to  improve  on  syllable  sorting,  a  program  must 
pull  out  speech  sounds  that  are  more  than  10  dB  below  the  competing 
sounds  from  another  speaker. 

5.2.  SORT  BY  MISTUNING  ALGORITHM 


A  second  special  case  which  allows  assigning  sounds  to  speakers  can 
occur  when  the  voice  signals  are  transmitted  by  single  sideband  (SSB) 
radio.  For  SSB  reception,  it  may  happen  that  one  voice  is  received  with 
an  amount  of  mistuning  differing  from  that  of  the  other.  The  use  of 
complex  correlation  then  gives  evidence  related  to  voice  frequency  shift, 
and  hence  to  voice  mistuning.  If  two  voices  are  being  received  with 
suitably  different  mistunings,  and  if  the  amounts  of  the  two  mistunings 
are  known,  then  the  voices  can  be  sorted.  Section  6.2  describes  a  program 
developed  during  the  present  project  which  automatically  estimates  pitch 
and  mistuning  for  single  voice  speech.  The  program  also  produces  histo¬ 
grams  in  graphic  form  that  can  help  characterize  noisy  and  mixed  voice 
signals . 


As  shown  in  Section  3.3,  the  complex  correlation  of  voiced  speech 
data  may  be  used  to  determine  what  raised  sinusoid  comb  best  fits  the 
speech  power  spectrum.  For  C(n)  the  complex  correlation,  and  T  the  value 
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of  n  (within  the  range  of  possible  pitch  periods)  that  maximizes  C(n), 
then  the  best  comb  F  is 


F(k)  =  0.5  +  0.5  eos(27TTk  +  d) 

where 

d  =  arctan(-Im  C(T)/Re  C(T))  . 


For  properly  tuned  speech,  d  should  be  zero.  This  is  because  the  peaks 
in  the  spectrum  are  all  harmonics  of  the  pitch.  The  comb  which  best  fits 
the  spectrum  will  therefore  have  a  maximum  (a  tooth)  at  zero  frequency. 
Deviation  of  the  lowest  frequency  tooth  from  zero  frequency  indicates 
speech  mistuning. 

A  pitch  period  of  T  sample  intervals  is  a  pitch  period  of  T/F  seconds 

s 

where  F  is  the  sampling  frequency  in  Hertz.  The  pitch  is  therefore  F  /T 
s  s 

Hertz.  Speech  mistuned  by  H  Hertz  is  therefore  mistuned  by  H/pitch  = 

HT/F  times  the  pitch.  This  is  2nHT/F  radians  in  the  cosine  term  of 
s  s 

F(k).  The  comb  itself  is  displaced  by  -d  radians.  Therefore  the  phase 
angle  between  a  comb  fitted  to  a  mistuning  of  H  Hertz  and  the  actual  comb 
is  (d  +  2ttHT/Fs)  modulo  2v.  When  this  phase  angle  is  small  (or  close  to 
a  multiple  of  2tt),  the  observed  data  is  consistent  with  the  hypothesis 
that  a  frequency  displacement  of  H  Hertz  has  taken  place.  This  is  the 
basis  of  the  sort  by  mistuning  algorithm. 


The  sort  by  mistuning  program  is  given  two  mistunings  to  check  for. 
It  then  goes  through  the  data,  and  for  each  time  window  determines  the 
bes'  comb  filter  fit.  It  checks  the  phase  of  each  comb  to  see  which 
mistuning  better  predicts  the  phase  angle  of  the  comb.  For  --.ach  of  the 
two  mistunings  there  is  an  output  channel.  Lach  data  interval  is  sent  to 
the  output  channel  whose  mistuning  better  predicts  the  comb  phase.  If 
neither  mistuning  predicts  the  comb  phase  well,  then  the  data  is  sent  to 
both  output  channels.  If  an  output  channel  does  not  receive  data,  then 
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the  gap  is  filled  with  white  noise.  The  sort  by  mistuning  had  the  usual 
parameters,  namely  6400  Hertz  sampling  rate,  512  point  (80  milliseconds) 
raised  cosine  time  windows,  and  time  window  overlap  of  50  percent. 

As  the  preceding  discussion  implies,  mistunings  can  best  be  used  to 
separate  data  if  they  predict  distinctly  different  comb  phase  angles,  most 
of  the  time.  This  is  most  likely  if  both  speakers  speak  with  about  the 
same  pitch,  and  the  difference  in  mistunings  is  about  half  of  the  common 
pitch.  Listening  tests  confirmed  that  under  this  condition  the  sort  by 
mistuning  algorithm  separates  voices  well. 

The  data  sorting  algorithms  separated  voices  well  under  their 
enabling  conditions,  but  there  was  clearly  room  for  improvement.  The 
white  noi' e  introduced  was  distracting,  and  occasional  absences  of  data 
were  noticeable.  The  idea  suggested  itself  to  apply  the  comb  filtering 
described  in  Section  3.  The  data  sorting  criterion  would  then  be  used  to 
switch  between  comb  (enhancement)  and  rake  (suppression)  filtering, 
instead  of  switching  between  data  and  noise.  Section  6  discusses  both 
this  algorithm  and  a  program  to  automatically  estimate  mistuning. 
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SECTION  6 

SYLLABLE  COMBING  AND  MISTUNING  ESTIMATION 


The  programs  described  in  Sections  3  and  5  were  combined.  The  result 
was  the  best  voice  processing  developed  during  the  project.  This  voice 
processing  was  used  tc  make  the  project  demonstration  tape.  The  voice 
processing  was  aided  by  an  automated  method  to  estimate  pitch  and  SSB 
mistuning. 


6.1.  SYLLABLE  COMBING  VOICE  PROCESSOR  AND  RESULTS 

Two  syllable-combing  programs  were  written.  Each  operates  on  data 
sampled  at  6400  Hz.  Each  uses  512  point  (80  millisecond)  raised  cosine 
(Hanning)  windows  overlapped  50  percent.  Each  produces  only  one  output 
channel.  Each  combs  (enhances)  or  rakes  (suppresses)  each  data  windov; 
depending  on  a  selection  criterion.  These  programs  are  named  FRCOMB  and 
MTCOMB,  and  are  listed  in  Appendix  B.  FRCOMB  selects  by  pitch  frequency. 
The  operator  tells  the  program  the  pitch  range  to  check  for.  Windows 
with  pitches  in  this  range  are  combed.  All  other  windows  are  raked. 

MTCOMB  selects  by  mistuning.  The  operator  tells  the  program  what  fre¬ 
quency  shift  to  check  for.  The  operator  also  tells  the  program  what 
percent  of  combing  to  use.  This  percent  combing  indicates  what  complex 
correlation  phase  angles  are  to  trigger  data  combing.  For  example,  for 
50  percent  combing  the  program  would  check  each  window  for  a  comb  phase 
angle  within  50  percent  of  tt  radians  (the  maximum  possible  angular  dis¬ 
tance)  of  the  predicted  phase  angle.  If  the  phase  angle  is  smaller  cKan 
the  limit,  then  the  window  is  combed.  Otherwise  the  window  is  raked. 

The  program  which  merges  two  voice  data  channels  is  called  WTMERG  and 
is  listed  in  Appendix  B.  A  second  voice  merging  program  was  written  which 
adds  white  Gaussian  noise  at  a  specified  gain  level.  This  program  is 
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called  WTNMRG,  and  is  also  listed  in  Appendix  B.  Listening  tests 
showed  that  the  voice  combing  algorithms  are  robust  against  white 
Gaussian  noise. 

The  syllable  combing  programs  separated  voices  well  when  their 
selection  criteria  were  valid.  Each  raked  (suppressed)  voice  was  still 
slightly  audible,  but  sounded  faint  and  whispery.  Occasionally  the  pro¬ 
gram  would  suppress  portions  of  the  voice  which  should  have  been  enhanced, 
but  the  traces  that  remained  helped  the  ear  fill  in  the  gap  and  maintain 
continuity. 

When  trying  to  enhance  received  single  sideband  data  it  helps  to 
know  the  mistuning.  Similarly,  when  trying  to  sort  by  pitch  it  helps  to 
know  speaker  pitch.  A  program  to  determine  pitch  and  mistuning  is 
discussed  in  the  next  subsection. 

6.2.  AUTOMATIC  SSB  MISTUNING  ESTIMATOR 

As  described  in  Section  5.2,  a  given  best-fit  comb  phase  is  con¬ 
sistent  with  many  possible  mistunings.  For  F^  a  possible  mistuning 
frequency  and  F^  the  pitch  frequency,  then  F^  +_  F^  is  also  a  possible 
mistuning  frequency.  Thus  given  one  data  window  there  are  multiple 
solutions  to  the  mistuning  estimation  problem.  However,  speaker  pitch 
will  generally  vary  over  time.  For  a  changed  pitch,  the  new  multiple 
solutions  will  tend  not  to  line  up  with  the  old  multiple  solutions, 
except  at  the  true  mistuning.  This  can  be  used  to  estimate  mistuning. 

Program  PMCES  is  a  pitch  and  mistuning  estimator,  and  is  listed  in 
Appendix  B.  The  program  operates  with  the  usual  parameters,  namely  a 
6400  Hertz  sampling  rate,  and  512  point  (90  millisecond)  raised  cosine 
(Hanning)  data  windows.  The  windows  are  overlapped  by  50  percent.  A 
program  was  also  tried  with  256  point  (40  millisecond)  windows.  No 
major  difference  in  program  results  was  observed. 
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The  program  compiled  a  two  dimensional  histogram,  representing 
estimated  mistuning  versus  estimated  pitch.  There  were  18  pitch  bins, 
representing  80  to  90  Hertz,  90  to  100  Hertz,  etc.,  through  250  to  260 
Hertz.  There  were  161  mistuning  bins,  representing  -402.5  to  -397.5, 

-397.5  to  -392.5,  etc.,  through  397.5  to  402.5  Hertz.  This  gave  an  18  x 
161  =  2898  bin  histogram.  For  each  time  window,  the  program  computed  the 
complex  correlation.  It  then  found  the  correlation  magnitude  peak  in  the 
range  of  80  to  256  Hertz.  From  the  peak,  the  program  estimated  speaker 
pitch,  and  also  estimated  the  multiple  solutions  to  the  mistuning  estima¬ 
tion  problem.  Then  the  program  incremented  the  corresponding  bins  in  the 
pitch-mistuning  histogram. 

The  histogram  increment  used  was  the  magnitude  squared  of  the  complex 
correlation  peak.  Using  this  increment  gave  relatively  less  weight  to 
windows  whose  spectra  did  not  have  the  comb-like  structure  typical  of 
voiced  speech.  This  increment  also  is  related  to  the  power  in  a  window, 
because  the  complex  correlation  is  proportional  to  speech  amplitude.  A 
program  was  also  written  which  incremented  the  histogram  by  window  power. 
Incrementing  by  the  correlation  peak  was  found  to  give  a  cleaner  plot 
when  applied  to  very  noisy  data. 

In  addition  to  incrementing  the  histogram,  the  program  also  calculates 
the  weighted  average  of  the  estimated  pitch  and  the  weighted  average  of 
the  square  of  the  estimated  pitch.  The  weighting  used  is  the  same  as  the 
histogram  increment.  When  all  the  data  is  processed,  the  program  computes 
and  prints  the  mean  and  standard  deviation  of  estimated  speaker  pitch. 

It  also  estimates  mistuning  by  finding  the  peak  of  a  one-dimensional  mis¬ 
tuning  histogram.  Finally,  it  plots  the  pitch-mistuning  histogram. 

Figures  6-1  and  6-2  show  some  plots  that  this  program  produced. 

Figure  6-1  results  from  processing  37  seconds  of  microphone  speech  of  a 
male  speaker.  Mistuning  is  plotted  horizontally  and  pitch  is  plotted 
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Figure  6-1  Microphone  Speech  Pitch  Versus  Mistuning 


vertically.  The  pitch  is  seen  to  vary  mainly  between  100  to  150  Hertz. 

On  each  line,  multiple  solutions  for  the  mistuning  are  visible.  However, 
the  only  mistuning  indication  common  to  all  lines  is  at  the  vertical 
axis,  which  represents  zero  Hertz  mistuning.  It  is  typical  of  these 
pitch-mistuning  plots  that  the  peaks  form  slanting  lines  for  false  mis¬ 
tuning  solutions  and  a  vertical  line  for  the  true  solution.  Sometimes 
this  property  can  be  used  to  determine  the  true  mistuning  by  eye,  for  very 
noisy  data  that  fools  the  automatic  mistuning  calculation. 

Figure  6-2  results  from  reception  of  amateur  single  sideband  (SSB) 
under  moderately  noisy  conditions.  It  is  for  ten  seconds  of  data  during 
which  the  speaker  said  "...  uniform,  W2PAU,  five  nine  zero  five,  five 
nine  zero  fi'^e."  This  plot  :.s  not  as  clean  as  Figure  6-1,  but  the  speaker's 
pitch  and  mistuning  are  clearly  visible.  The  speaker's  pitch  was  110  to 
190  Hertz.  The  frequency  shift  of  the  received  signal  was  145  Hertz, 
indicated  by  the  vertical  line  of  peaks. 

These  results  show  that  sorting  by  pitch  or  mistuning  does  not  have 
to  be  based  on  guesswork.  An  automated  method  is  available  for  estimating 
pitch  and  mistuning. 


SECTION  7 

CONCLUSIONS  AND  RECOMMENDATIONS 


As  Sections  2  through  6  have  shown,  the  main  positive  results  of  the 
present  effort  were  based  on  the  use  of  complex  correlation.  Adaptive 
lattice  filtering  did  not  lead  to  useful  results.  The  complex  correlation 
allowed  automatic  adaptive  design  of  comb  filters  to  enhance  and  suppress 
data.  The  complex  correlation  also  gave  estimates  of  speaker  pitch,  and 
could  be  used  to  estimate  speaker  frequency  shift.  It  also  turned  cut  to 
be  robust  against  noise  and  distortion. 

The  voice  separation  algorithms  developed  during  the  project  apply  to 
special  cases.  The  main  obstacle  to  treating  the  general  case  was  the 
lack  of  a  general  method  to  identify  to  which  speaker  the  various  sounds 
in  mixed  speech  belong.  Identifying  by  adaptive  lattice  filtering  coeffi¬ 
cients  (Section  2)  was  found  to  be  impractical.  Identifying  by  power 
level  (Sections  3  and  4)  was  also  found  to  be  unsatisfactory.  Voice 
separation  was  achieved  for  the  special  cases  of  distinct  speaker  pitches 
and  distinct  frequency  shifts  (Sections  5  and  6). 

A  program  for  .automatically  estimating  speech  frequency  shift  was 
developed  in  this  project.  This  may  be  a  new  result.  In  itself,  it  may 
aid  in  improving  intelligibility,  or  at  least  it  may  be  used  to  improve 
single  sideband  reception  quality.  The  estimator  can  be  used  to  estimate 
the  mistuning  of  baseband  signals,  and  thereby  direct  their  retuning. 

Whether  or  not  the  techniques  developed  will  aid  in  any  particular 
radio  reception  problem  may  depend  on  the  details  and  specifics  of  the 
problem.  There  are  thus  two  ways  to  go  in  the  future.  One  way  is  to 
take  the  algorithms  to  the  problems.  By  this  is  meant  developing  software 
(and  possibly  host  hardware)  which  is  convenient  for  data  analysts  to  use 
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and  apply  to  their  particular  problems.  The  algorithm  developed  to 
separate  speech  can  be  programmed  to  run  in  real  time,  if  a  dedicated 
FFT  or  array  processor  is  used. 

A  second  way  to  proceed  in  the  future  is  to  bring  the  problems  to 
the  algorithms.  By  this  is  meant  supplying  some  target  data  for 
processing  with  the  existing  software,  perhaps  with  modifications  to  better 
fit  the  data.  Each  approach  has  its  own  advantages.  Experienced  data 
analysts  are  more  likely  to  appreciate  the  results  of  signal  processing. 

On  the  other  hand,  the  second  approach  is  much  less  costly,  and  would 
allow  signal  processing  by  personnel  familiar  with  the  algorithms  developed. 
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APPENDIX  A 

SSB  MISTUNING  SIMULATOR 


A  program  was  written  to  perform  digitally  the  frequency  shifting 
needed  to  simulate  single  sideband  (SSB)  mistuning.  The  program  was  also 
useful  for  correcting  the  mistunings  of  actual  SSB  radio  data.  The  pro¬ 
gram  is  named  WMTSSB,  and  is  listed  in  Appendix  B. 

The  heart  of  the  program  is  a  90  degree  phase  shifter,  also  known  as 
a  Hilbert  transformer.  This  is  a  finite  impulse  response  (FIR)  digital 
filter.  It  was  designed  using  an  FIR  filter  design  program  (using  the  Remez 
exchange  algorithm)  which  was  taken  from  Rabiner  and  Geld  [4].  A  real 
waveform  is  always  made  up  of  positive  and  negative  frequency  components. 
However,  an  imaginary  part  can  be  attached  to  a  pure  real  waveform,  in 
order  to  make  the  waveform  have  only  positive  frequency  components.  This 
process  simulates  in  baseband  the  sideband  suppression  filtering  normally 
done  at  IF  in  an  SSB  transmitter. 


Lower  sideband  suppression  may  be  done  by  applying  a  linear  filter 
with  frequency  response  G(e~’u))  =  0.5  +  0.5jH(e')tJj)  where 


H(eja)) 


0  <  a)  <  v 


TT  <  0)  <  2  TT 


Here  we  are 
transform  of  the 
Gold  show  [4,  p. 


working  with  a  sampled-time  system,  so  H(e~*w)  is  the  z 
impulse  response  h(n)  with  z  =  e^u.  Then  as  Rabiner  and 
71] 


h(n) 


2  sin2(Trn/2) 
im 


n  i  0 

n  =  0 
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This  is  the  ideal  digital  Hilbert  transform  filter.  Rabiner  and 
Gold  include  in  their  book  (pp.  194  to  204)  a  program  that  can  design 
finite  impulse  response  (FIR)  approximations  to  the  ideal  Hilbert  trans¬ 
former.  This  filter  design  program  was  implemented  by  PAR  on  a  previous 
project . 

The  filter  design  program  was  used  to  design  an  equiripple  approxi¬ 
mation  to  Hie^),  for  the  passband  of  frequencies  from  0.03  to  0.47 

F  ,  where  F  is  the  sampling  frequency  (6400  Hertz  for  the  present 
s  s 

project).  This  is  a  passband  of  192  to  3008  Hz.  A  31  tap  filter  was 
designed  which  turned  out  to  be  wj thin  +_  0.21  dB  of  the  ideal  within  the 
passband.  All  the  even  numbered  taps  were  zero.  The  odd  numbered  taps 
are  given  in  Table  A-l.  The  taps  of  the  table  are  related  to  the  ideal 
h(n)  by  g(n+16)  approximates  -h(n). 


Table  A-l 

Hilbert  Tran 

sformer 

Tap  Weight 

Tap 

Weight 

Tap 

g(l) 

0.02050 

= 

-g(31) 

g(  3 ) 

0.02134 

= 

-g(29) 

g(5) 

0.03265 

- 

-g( 27 ) 

g(7) 

0.04876 

= 

-g( 25 ) 

g(9) 

0.07296 

= 

-g( 23 ) 

g(ll) 

0.11398 

= 

~g( 21 ) 

g(13> 

0.20402 

= 

-g(19) 

g(15) 

0.63385 

= 

-g(17) 

The  frequency  shifting  program  used  the  Hilbert  transformer  to  zero 
the  negative  frequency  components  of  the  signal.  Next  it  doubled  the 
sampling  rate  for  the  (now  complex)  signal.  This  was  done  so  that 
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frequency  shifting  would  not  result  in  aliasing.  For  example,  a  sinusoid 
at  3  kHz  shifted  up  400  Hertz  to  3,4  kHz  would  appear  not  to  be  shifted 
at  all,  if  the  sampling  rate  were  only  6.4  kHz.  This  is  the  well  known 
"aliasing"  effect  for  sainpled-time  signals.  The  sampling  rate  was  doubled 
by  first  inserting  a  zero  between  each  pair  of  consecutive  samples.  Then 
the  signal  was  passed  through  a  six  pole  Butterworth  lowpass  filter, 
which  smoothed  the  data.  The  filter  cutoff  was  0.2  of  the  new  (doubled) 
sampling  rate.  This  is  2660  Hertz.  Next  the  signal  was  multiplied  by  a 
complex  exponential,  producing  a  frequency  shifted  signal.  Only  the  real 
part  of  the  result  was  saved.  Next  th^.  signal  was  lowpassed  with  the  same 
Butterworth  filter  used  previously.  This  prevented  aliasing.  Finally, 
the  sampling  rate  was  halved,  that:  is,  returned  to  its  original  value. 


To  summarize ,  the  steps,  in  the  frequency  shifter  were 

1.  Attach  j  times  the  Hilbert  transform  to  the  signal. 

2.  Double  the  sampling  rate  by  interpolating  zeroes. 

3.  Smooth  the  data  by  lowpass  filtering. 

4.  Multiply  by  a  complex  exponential  to  shift  frequency. 

5.  Discard  imaginary  part. 

6.  Lowpass  again  co  prevent  aliasing. 

7.  Discard  alternate  samples. 

The  Hilbert  transformer  was  implemented  so  that  tap  16  represented 
time  zero  of  the  impulse ,  response .  This  made  it  a  "non-realizable."  filter, 
which  predicts  the  future.  This  was  done  to  avoid  introducing  delay  in 
the  <5 1 'jti a  1  uhon  it  wi  +-+  the  imapinarv  sienal  that  was  attached 

to  it . 


Figure  A-l  shows  in  graphic  form  the  result  of  frequency  shifting  by 
+200  Hertz.  Lines  in  the  figure  are  organized  in  groups  of  three.  Each 
line  consists  of  four  frames,  with  each  frame  of  each  line  normalized. 
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Each  frame  has  256  points  (40  milliseconds)  of  data.  The  first  line  of 
each  group  shows  the  original  data.  The  second  line  of  each  group  shows 
the  imaginary  part  attached  to  the  data.  Note  that  this  is  in  phase 
quadrature  with  the  original  data.  That  is,  the  second  line  of  each 
group  has  a  zero  crossing  where  the  first  line  has  a  peak,  and  vice  versa. 
The  third  line  shows  the  result  of  shifting  the  data  up  200  Hertz.  Plots 
such  as  this  nelped  verify  the  correct  operation  of  the  program. 


APPENDIX  B 
COMPUTER  SOFTWARE 


The  software  to  be  presented  was  developed  to  run  on  a  Digital 
Equipment  Corporation  PDP-11/45  computer  with  a  Tektronix  4014-1  alpha- 
numerics/graphics  terminal  and  a  Tektronix  4631  hard  copy  unit.  The 
software  runs  under  the  RSX-11M  operating  system,  version  3.1.  All  the 
software  developed  during  the  project  was  compiled  under  Fortran  IV 
Plus.  The  software  is  believed  to  be  compatible  with  Fortran  IV,  but 
this  has  not  been  checked.  Several  subroutines  were  used  for  input  and 
output  that  were  not  developed  during  the  project,  and  are  not  listed. 
These  are  READF,  WRITEF,  PLOTS,  PLOT,  and  LINE. 


READF  and  WRITEF  are  fast  programs  to  read  and  write  data  to  and 
from  disk.  Their  calling  formats  are 


CALL  WRITEFCLUN,  BLKNO ,  BUF,  SIZE,  IOSB) 
CALL  READF  (LUN,  BLKNO,  BUF,  SIZE,  IOSB) 


where  LUN 

BLKNO 

BUF 

SIZE 

IOSB 


is  the  logical  unit  number  (integer  *2) 

is  the  virtual  block  number  (integer  *4) 

is  the  data  buffer  (vector  cf  integer  *2 ) 

is  the  data  clock  size  in  bytes  (integer  *2) 

is  the  I/O  status  block  (2  element  integer  *2  vector) 


PLOTS  is  called  with  no  arguments.  It  erases  the  Tektronix  4014 
screen,  puts  it  in  graphics  mode,  and  sets  the  plotting  origin  at  the 
lower  left  corner.  PLOT  is  used  in  the  project  software  to  move  the 
plotting  origin.  The  calling  format  is 


CALL  PLOT  (X,  Y,  -3) 


where  X  is  the  horizontal  shift  of  the  origin  (type  real) 

Y  is  the  vertical  shift  of  the  origin  (type  real) 

LINE  is  used  to  plot  a  line  of  data.  That  is,  LINE  plots  a  series  of 
points  and  connects  them  with  straight  lines.  The  calling  format  is 

CALL  LINE  (XBUF,  YBUF,  NETS,  1,  0,  0) 

where  XBUF  is  the  vector  of  X  coordinates  (vector  of  reals) 

YBUF  is  the  vector  of  Y  coordinates  (vector  of  reals) 

NPTS  is  the  number  of  points  to  plot  (integer  *2) 

XBUF  and  YBUF  must  each  have  NPTS+2  components.  XBUF  (NFTS+  i  )  and  YBUF 
(NPTS+1)  each  contains  the  starting  value  of  the  line  (usually  the  minimum 
value  in  the  line).  XBUF  (NPTS+2)  and  YBUF( NPTS+2)  tell  the  number  of 
coordinate-units  per  distance-unit. 

The  format  of  the  project  waveform  files  was  as  follows.  Each  file 
consisted  of  a  series  of  512-byte  blocks.  The  first  block  of  each  file 
was  a  header  block.  This  was  followed  by  a  variable  number  of  data  blocks, 
each  containing  256  integer  *2  words.  The  header  block  format  in  integer 
*2  words  was  as  follows. 


Word( s ) 

Value 

Mean ing 

1 

1 

One  channel  of  data 

(2,3) 

(6400,0) 

Samj ling  rate  6400  Hertz 

(4,5,6) 

- 

ttart  time  in  hours,  minutes, 
seconds 

(7,8,9) 

- 

Stop  time  in  hours,  minutes, 
seconds 

(10-13) 

- 

Not  used 

14 

0 

Data  is  integer  *2 

(15,16) 

- 

Not  used 

(17,18) 

(NL0,NHI ) 

There  are  65536*NHI+NL0  data 
blocks. 
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Not  all  programs  maintained  the  start  and  stop  times.  Words  17  and  18 
sometimes  appear  in  the  programs  as  word  9  of  an  integer  *4  array. 

The  software  produced  during  th°  program  was  written  using  structured 
programming,  with  statement  indenting  to  show  program  structure.  Use  of 
this  method  was  found  to  greatly  aid  program  readability  and  reliability. 

The  program  listings  have  hand-drawn  f low-of-control  arrows  added. 
This  makes  each  listing,  in  effect,  its  own  flowchart.  Programs  are 
listed  in  the  alphabetical  order  of  the  main  programs  in  which  they 
appear.  The  exceptions  are  files  CFFT  and  HILB.  File  CFFT  contains  the 
complex  fast  Fourier  transform  (CFFT)  subroutine.  File  HILB  contains  two 
subroutines  used  by  the  WMTSSB  program.  These  are  HILB  and  LWPS,  the 
Hilbert  transformer  (90  degree  phase  shifter)  and  Butterworth  lowpass 
filter  respectively. 

Table  B-l  gives  a  table  of  contents  for  the  listings.  SBR  indicates 
that  a  program  is  a  subroutine. 


s 


Table  B-l  Guide  to  Software  Listings 


Program 


Function 


CCORPL 

CFFT 

FRCOMB 
HI  LB 

MTCOMB 

PMCES 

VCRFSP 

WMTSSB 

WTMERG 

WTNMRG 


Complex  correlation  and  plot 

Complex  fast  Fourier  Transform 
( FFT )  SBR 

Comb  by  pitch  frequency 

Hilbert  transformer  SBR,  lowpass 
SBR 

Comb  by  (SSB)  mistuning 

Pitch-mistuning  estimator 

Voice  lattice  filter,  threshold, 
reconstruct 

"Write  mistuned  SSB"  (Frequency 
shifter) 

Does  weighted  merge  of  waveform 
files 

Does  weighted  merge,  adds  white 
Gaussian  noise 
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c — ccorpl  by  bob  dick 

C - FROM  TNSPPl  14  AUG  SO 

C - COrtPLKX  CORRELATION  MAGNITUDE  PLOT 

C - ONE  SECOND  OF  BATA  PER  PA6E. 

BYTE  FNAM(30)»ANSW 

DIMENSION  DATIN( 2541 iWNDM<256>»TR(25N))TJ (256) 
DIMENSION  AMPL<5?MXCi52MDUF<25fti 
INTEGER  14  IRECrlOSBtlBFdi 
EQUIVALENCE  (IBF.IBUf) 

DATA  DATIN/256W./PI/3.14l5WaS/ 

— >TYPE  10010 

ACCEPT  10020)1  NTHt(FNAM(RH/  tMi-'iil  NTH/ 
FNAH(LNTH+t )r0 

OOP£N(UNIT3)  »NAME*FNANiTYPE='OlD'  > 

1  BUFFERCOUNT =- 1 i ERR-9000 1 RE  ADONl i / 


IRECM 

CALL  READF(l»IREi  iIW'F»512)I0SB> 
NlBLKS=IDEt9) 

TYPE  10030f N1BLKS 
ACCEPT  10040)NFBLK 

IFTNFBLKd  E. 0)G0  TO  9000 - >  AB6RT 

ACCEPT  10040. NMBLRS 
NLBLK-NFBLK+NMBLKS-J 

IF(NIBLK,LT.NFBLK)G0  TO  9000 — ►  ABORT 
1F(NLBLK,GT.N1BLKS/G0  TO  9000— >  ABORT 
C - SETUP  HANNING  BATA  WINDOW 

flO  IDXS1*256 

WNDW(1DX)=0. 5-0, 5IC0S, Pill jx,  (26./ 

10  CONTINUE 

C - SET  X  COORDS  FOR  AMPL  LINE 

DO  20  IDX*t)50 
f  XC<IDX)*5MDX 
20  CONTINUE 
XC(51)=1 
XC<52)=1 */0. 23 

C - SET  1  CHORDS  FOR  CQRR  LINES 

CALL  STYC 
NFIRST*NFBLKI2-1 

C . SNIP  FIRST  HALF-WINDOW 

CALL  BA128(DATIN)NF1RS,) 

NMIN=NF1RST)1 

LAST-NLBI.K42 

C - DO  UNTIL  (NO  MORE  DATA)  PLOT  PAGES 

30  CONTINUE 

C - A--DONE  ONCE  PER  PAGE 

j  CALI.  PLOTS 
DO  40  IDX-1.50 
f  AMPLUDX)=0. 

40  CONTINUE 

DO  90  IxLN=l >50 

C . A--DONE  ONCE  PER  CORRELATION  LINE 

C . . GET  NEW  DATA 

CALL  0A128(DATIN)NHIN) 

PW-O. 

DO  50  IDX= 1 » 256 

ji  A  TR ( IDX ) =DAT IN ( IDX ) 4WNDW l I DX ) 
N3PWFTR(IDXJIDATIN<  idx  ■ 
TI(1DX)*0. 
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50  CONTINUE 

AHPL<IXLN)=S»T(PU) 

CALL  CFFTU.TR/TI.B) 

: - FORM  AMPLITUDE  SPECTRUM 

DO  60  IDXd.128 

f  TR(IDX)  *SVRT  { TR(  IDX )  *TR(  IDX )  HI  <  IDX )  III  ( IDX ) ) 
|  TI(IDX)-0. 

60  CONTINUE 

C - ZERO  NEGATIVE  FREQUENCY  PORTION 

DO  70  IDX=129,256 
f  TR(IDX)=0, 

|  TI(IDX)"0- 
70  CONTINUE 

C - — TAKE  INVERSE  FFT 

Cfcl  CFFT,-1/TR»TI.8) 

C - FORft  MAGNITUDE  OF  COMPLEX  CORRELATION 

DO  80  IDX=1.128 

A  A  ^  TRTIDX)=SORT(TR(IDX)*TR(IDX)fTJ{IDX)tTI(IDX)) 
80  CONTINUE 

C - PLOT  A  CORRELATION  LINE 

CALL  CCRPLT(TRiIXLN) 

NMINcNHINtl 

C . — IFTNQ  MORE  DATA)  EXIT  CORR  LINE  LOOP 

IF(NNIN.GT.LAST)GO  TO  100 
TO  CONTINUE 

100  CONTINUE 

C - PLOT  ANPLITUHE  LINE 

CALL  FCTRIAMPL.52) 

AMPL ( 52 )=AMPL (52)70.8 
Y0RG=9.9 

CAM  PLOKO.iYORGf-j) 

CALL  UN£(XC.AMPI.>  50. 1.0.0) 

C - -|— REPEAT  1 INE  TO  FILL  BUFFER 

I  CALL  l INE<XC« AWL. 50. 1.0.0/ 

A  yomh-yorg 

CALL  PLOT(O..YORG»-3) 

C - HAD  FOR  SIGNAL 

ACCEPT  10050. ANSU 

C - IF(X  TYPED1  EXIT  PAGE  LOOP 

IF<AN3H.E0.'X‘.GQ  TJ  1)0 

C - IF<MORE  DATA)RE,;>frtT  PAGE  LOOP 

IF<i<MIN  LE.LASDGO  TO  30 
110  CONTI,' NJF 

9000  CQN!!MJE< - /ABORT 

STOP 

10010  FORNATC  COMPLEX  CORRELATION  PLOTTER  INPUT  FILE?') 
10020  F0RMATIQ.30A1) 

10030  F0RMATUX.I8. '  B10CKS.  GIVE  FIRST  <CR>  HON  MANY') 
10040  FQRMAT(IB) 

10050  FORMAT(Al) 

_ END _ 

SUBROUTINE  BA128(DAT.KALL) 

C - IT  DIC*  8  AUG  80 

C - ADVANCES  DATA  128  PTS/CALL 

DIMENSION  DATi256).IBUF<256) 

INTE8ERI4  IREC.IOSF 

— >  IF(.N0T.(XALL.Fa.2«KALL/2)))6O  TO  20 

C - j— HERE  CALL  IS  EVEN.  USE  CURRENT  BUFFER 
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PACE  3 


uu 

t 


DO  10  IDXM»128 

DAT ( IDX )=DAT  < IDX+1 28 ) 
DAT(IDX+128)*IBUF(IDX+128) 
10Y  CONTINUE 
GO  TO  40 
20  CONTINUE 

C - 1 — HERE  CALL  IS  ODD,  USE  NEW  BUFFER 

IREC-KAi.L/2+2 

CALL  R£ADF(MRECrIBUF.512fIQSB) 
DO  30  IDX=l f 128 

DAT<IDX)CDAT(JPXF12S) 

,  CAT  <  IDX+128) =1  BLfF  ( IDX ) 

30  Y  CONTINUE 
40  CONTINUE 
—  RETURN 
END 


UL 

t 


SUBROUTINE  CCRPLT(DATAiLlN) 

C - BY  BOB  DICK 

C - FRON  SPCF1T  14  AUG  80 

C - PLOTS  A  LINE  OF  COMPLEX  CORRELATION 

DIMENSION  DATA(130)»YC(130) 

— >  CALL  FCTR( DATA. 130) 

DATA(130)=DATA(130)/O.A9 

X0RG=(50-LIN)*0.21 

C - 0.23449+0.69-11.96  INCHES  HORIZ 

CALL  PLOT  ( XORG » 0 » » --3 ) 

CALL  LINE(DATA.YC.128il»0<0) 
X0R6=-X0RG 


CALL  PL0T(X0RG,0.,-3) 

RETURN 
ENTRY  STYC 

->D0  1010  IDX-lsl28 
t  YC(JI'X)=128-IBX 
1010  CONTINUE 
YC(129)=0. 

C - 128  POINTS  AT  13/INCH  IS  9.8  INCHES 


YC<130)=13. 

<-  RETURN 

END _ 

SUBROUTINE  FCTR(FNCrlfUM) 

DIMENSION  FNC(2S8) 

->  IF( .NOT.  <IDIM.GE.4>)G0  TO  900 - 

IF(.NOT.(IDIM.LE.258))GO  TO  900 - 

FMN=FNC(1) 

FMX=FNC(l) 

DO  10  IDX=2» IDIM-? 

IF ( FHN , 6T , FNC ( I PX ) ;FNN=FNC  < IDX ) 

I  IF ( FMX . L  T . FNC ( IDX ) ) FMX=FNC  t IDX  > 

10  CONTINUE 

FNC(IDIM-1)=FMN 

FNC(IDIM)=FMX*FMN 

900  CONTINUED - ABORT 

< —  RETURN 
END 


ABORT 

ABORT 
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SUBROUTINE  CFFT<MOO£,X*,XJ,N) 

C - »r  BOB  PICK,  ENTERED  14  f£B  SO, 

C - "FAST  FOURIER  TRANSFORM  WITH  COMPLEX  llftJT 

C - mKMORMh  HODf-*BACK«ARD 

rill  "SrilS  VECT0R’  5(1  IS  IHAC'  dimension  is  2*«f, 

r  IW  TIW£  ALG0RITHrt’  &KHHE1H  AND  SCHAFER 

C  FIG  6,10  P  297,  BUGGT  PROGRAM  ON  P  331. 

DIMENSION  XR(1024),XI(1024) 

BATA  Pt/3, 14159765358979/ 

->M=2tlM 

IW2*H/2 

C  ORDER  DATA  ST  SIT-REVERSED  JMDFX 
MMJ=N-1 
I»TRV=1 

M  40  IIiX=l,MMl 

If ( ,NOT i ( ID) ,LT, IBTRV) )G0  TO  JO 

wamutm) 

WMuimv} 

mmmmam 

Xi(IBTRV>=XI(IDX> 

xRum-im 
XKIDXJsTMPI 
10 1  CONTINUE 

"INCREMENT  rH£  BIT-REVERSED  INDEX 
IPTR=NV2 
GO  TO  30 
:0  I  CONTINUE 

A  IBTRV=IBTRV-IPTR 
f  IPTR=lPTR/2 
CONTINUE 

If (IBTRV, GT,IPTR)po  TO  20 
,  I8TRV=IBTfiVHPTR 
40  CONTINUE 

C - PERFORM  THE  BUTTERFLIES 

LEAP*! 

DO  70  ISTG=1 ,M 
JUMP*LEAP 
LEAP=L£AP« 

COI/NKW 
TWSTR'CO’.'IPI/COUNT) 
r«STl*SIW(PI/COUNT) 

If(MaPE.liE.O)TNSri--TySTl 
TURNR-1. 

TURN  ho, 

DO  60  IDX-h  JUMP 

DO  50  IPTR*JDX*NiLEAP 
Jpn»«jpT«.uttp 

TMPR=XR( JPTR; FTURNR-XI ( JPTRJITURNJ 
TMPi-AKl  JPTRIBTURNUxit  JPTRXIURNR 
XR(JPTR)=XRiIPTR)-TMPR 
XJ(JPTR)*XIMPTR)-TKPI 
.  XRUPTRhXRUPTRHTMPR 

j  XI<IPTR)-XI<IPTRHWl 
CONTINUE 

TMPR»TURNR*TNSTR'TURNJ*T«STI 
TURNhTURNRiTNSTITTURNilTWSTR 
T(JRNR»TNPR 
CONTINUE 


30  i 

t 


50 


H 


60 
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C - FRCOHB  BY  BOL  DICK 

C - LAST  UPDATE  10  SEP  80 

C - FROM  MTCOMB  4  SEP  80 

C - COM  FILTERS  VOICE  DEPEND INC  ON  PITCH  FREQUENCY 

C - DATA  ARRAYS! 

C - CCR.CCI!  UN  NAG  SPECTRUAt  THEN  COMPLEX  CORRELATION 

c — datin:  input  data  before  windowing 
c — ibuf:  input  buffer 

c - K8UF!  HEADER)  THEN  OUTPUT  BUFFER 

c — obuf:  data  overlap  output  buffer 

C - TRfTI!  WINDOWED  DATA)  THEN  F FT)  THEN  FILTERED  DATA 

C - WNDW:  RAISED  COSINE  (HANNINO)  WINDOW 

BYTE  F1NAH(30))F2NAM<30> 

DIMENSION  PATIN(512>)WNUy<512>)TR(512))TI<512> 
DIMENSION  CCRI512) )CCI (512) )0BUF(25A) 

DIMENSION  1BUF(256))KBUF(256) 

INTEGERI4  IREC)IOSB)IBF(?i)KBF(?) 

EQUIVALENCE  < I BF  r I BUF  >  > ( KBF , KBUF ) 

DATA  KBUF /256B0/DAT IN/51 210 . / 

DATA  OBUF/256IO . /RAKE ) COMB/ -0 .5)0.5/ 

DATA  PI.TWPI/3. 14159265)6. 28318531/ 

DATA  SAMPFQ/6400. / 

C - GET  INPUT  FILE)  NUMBER  OF  BLOCKS 

TYPE  10010 

ACCEPT  10020 )  LNTH  t  ( F1NAM  ( KH ) .  KH=  1 )  I.NTH) 
F1NAM(LNTHE1)=0 

OOPEN ( UNI T=1 ) NAME =F 1 NAM t T YPE= ' 01 D ' , 

1  BUFFERCOUNTs-1 ,ERR=9000) READONLY) 

IREC*1 

CALL  READF(1 )IREC) IBUF )512) IOSP) 

N1BLKS=1PF(9) 

TYFE  10030)N1BLKS 
ACCEPT  10040)NMBLKS 

IF(NMBLKS.LT ,2)60  TO  9000 - ►ABORT 

IF (NMBLKS.GT.NIBLKS)GO  TO  9000 - 9*  A  BOfcT 

C - -GET  OUTPUT  TILE)  WRITE  HEADER 

TYPE  10050 

ACCEPT  10020,LNTH)(F2NAM(KH))KH-1 >LNTH) 
F2NAM(LNTH+1)=0 

OOPFN(UNIT=2)NAME=F2NAM)TYPE='NEW'« 

1  BUFFERCOUNT * - 1  >  ERRS9000 ) 

KBUF<1)=1 
KBUF  (2)  =6400 
SEC=NMBLKS/25. 

MIN=NMBLKS/1500 
ISEC=SEC-60.4MIU 
KBUF(B)=MIN 
KBUF(9)=ISEC 
K8F(9)=NMBLKS 

CALL  WRITEF(2)IREC)KBUF)512)I0SB> 

C - SET  UP  DATA  WINDOW 

DO  10  JDX--1.512 

f  WNDW<IDX)*0.5-0.5*COS(PI*IDX/256.) 

10  CONTINUE 

C - 6ET  PITCH  RANGE  TO  COMB  FOR 

TYPE  10060 
ACCEPT  10070iPTl.OW 
ACCEPT  10070>P?HIGH 
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C - COMPUTE  PITCH  PERIOD  RANGE  FROM  PITCH  RANGE 

IF(,NOT.(PTLON.LT,PTH1GH))GO  TO  MOO - ^ ABORT 

PERHI*2S6. 

IF(PTL0W.6T,0.  >P£RHI«K400,/PTL0W 

IF(,H0T.(PTHI6H.GT.0.))G0  TO  9000 - > ABORT 

P£RL0*4400,/PTHI6H 

C - ZERO  OUTPUT  BUFFER 

DO  15  !DX=J,256 
f  K8UF(IDX)=0 
IS  CONTINUE 

C - READ  FIRST  DATA  RECORD  BEFORE  MAIN  LOOP. 

NNINS1 


C - + 


20 
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CALI  TISADVUIATIN'NfllN) 

C - LOOP  NUNBER.OFJtECORDS-l  TIMES. 

DO  80  NMIN=2»NMBLKS 

C— -^~HAIN  LOOP.  DONE  ONCE  PER  DATA  RECORD, 

■--TAKE  TINE  WINDOWS  WITH  50  PERCENT  OVERLAP. 

CALL  TFSADV(DATIN*NMIN) 

.DO  20  JDX*1>512 

TR<  JDX ) SDAT IN  < IDX ) IWNDW ( I DX ) 

T I ( I DX ) =0 . 

CONTINUE 

-FORM  AMPLITUDE  SPECTRUM 
CALL  CFFT (1»TR»TI»9) 

DO  30  IDX«1'2S6 

CCR(IDX)«SQRT<TR<IPX)ITR(IDXnTI<IDX)*TI(IDX)> 
CCI(IDX)=0. 

CONTINUE 

-SET  NEGATIVE  FREQUENCY  PORTION  TO  ZERO 
DO  40  IDX°257r512 
CCR(IDX>*0. 

CCI< IDX)=0. 

CONTINUE 

-TAKE  INVERSE  FFT 
CALL  CFFT(-l»CCR»CCIi?) 

-SEARCH  FOR  PEAK  FROM  I.OW  TO  HIGH  TIME  LIMIT 
PEAK=0 . 

IDPK*?5 

DO  SO  IDX*25»80 

TEMP»CCR ( IDX ) ICCR  < I DX ) +CC I ( IDX ) ICCI ( I DX ) 

IF(, NOT. (PEAK, IT, TEMP) )G0  TO  45 
I  PEAK-TEMP 
f  IOPK-IDX 
CONTINUE 
CONTINUE 

FIND  QUADRATIC  INTERPOLATION  PEAK, 

BSGsCCR( IDPK-1 )ICCR( IDPK-1 )ICCI ( IDPK-1 )BCCI ( 1DPK-1 ) 
USQ*CCR< IDPKFI )<CCR(IDPK+1 )FCCI (IDPK41 >ICCI ( IDPK41 ) 
X*l. 

IF(PEAK,OT,USQ)X-0. 

IF(BSQ, GE, PE AK)X*-1. 

IF<X,EQ,0.)X*0,5t(USQ-BSQ)/(2.#PFAK-BSQ-llSfl) 
TAU«IDPK-1FX 

INTERPOLATE  REAL  AND  IMAG  PARTS  (VIA  QUADRATIC) 
CF1M.5*(CCR(IDPKFI)-CCR<IDPK-1>) 
CF2*0,5t(CCR(IDPK-l )+CCR< IDPK+I ) )-CCR( IDPK) 
PARTR*CCR( IDPK ) M t CF1+XICF2 ) 
CFl*0,3t(CCI(IDPKFl)-CCI(IBPK-J)) 
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CF2=0.5*<CCI  ( WPX-mCCI<IBPKH))-CCI !  J.DPK) 
PttTI=CCI<lDP*HX*<CFl+X«CF2> 

PARTHrSflRT (PAR TRIP ARTREPARTItPART I ) 
IF<PARTN.LE,0,)PARTH=1. 

C . -COMB  OR  RAKE  ACCORDING  TO  PITCH 

SWITCH=CONB 

IF(TAU.IT,P£RLG)SWITCH=RmKE 
IFdAU.GT .PERHI )SWITCH=RAKE 
ONEG=PllTAU/256. 

DO  70  lDX=lt2S6 

A  CS= ( P ARTRICOS ( ONEG I 1 DX ) +PAR T II S IN ( ONEGI I DX ) ) /P ARTK 
FILT=0.5+SWITCHICS 
JDX=513-IDX 
TR(IDX)*FILT*TR(IDX> 

TI<IOX>=FILT*TI<IDX) 

TR(JDX)=FILTtTfi<JOX) 

Tl< JDX)=FILTITI< JDX) 

70  CONTINUE 

C - TtfE  INVERSE  FFT 

A  CALL  CFFT(-l»TR»TIi9> 

C - OVERLAP  RESULTING  TINE  FUNCTIONS 

NH0UT*NHIN-1 

CALL  LAPOUT(TR*OBUFf?»NNOUT) 

80  CONTINUE 

C . WRITE  FINAL  DATA  RECORD  AFTER  MAIN  LOUP, 

DO  90  IDX=1 » 256 

f  TR(inX)=0. 

90  CONTIMUC 

CALL  LAPOUT ( TR i OBUF 1 2 1 NNBLKS ) 

9000  CONTINUE < - * - ABOUT 

< —  STOP 

10010  FORHATt'  PITCH  FREQUENCY  COHBFR  INPUT  FILE?') 

10020  F0RNAT(0i30Al ) 

10030  FORHATUX* J8» '  BLOCKS.  HOW  MANY  DO  YOU  WANT?') 

10040  FORMAT ( 18) 

10050  FORMAT < '  WHAT  OUTPUT  FILE?') 

10060  FORNATt '  WHAT  PITCH  RANGE? (USE  .  AND  2  LINES)') 

10070  F0RNAT(F10.2) 

_ END _ 

SUBROUTINE  TFSADV(X.ICALL) 

C - BY  BOB  DICK.  REVISION  OF  14  MAR  80. 

C - ADVANCES  BUFFER  256  F0INT5  PER  CALL, 

DINENSION  X(512bIBUF!256> 

INTE6ERI4  IRECiIOSD 
IREC=ICALL+1 

CALL  READF <  1 .  IREC *IBUF.512.  IOSB ) 

DO  10  IDX-1.256 
f  XlIHX)*X(IDXE256) 

|  X<IDX4256)=IBUF(IDX) 

10  CONTINUE 
RETURN 

_ END _ 

SUBROUTINE  LAPOUT ( DAT .BUF (LUNrKALL) 

C - OUTPUTS  DATA  WITH  50  PERCENT  OVERLAP. 

C - INITIALLY  IUF  SHOULD  BE  ALL  ZERO. 

C - DOES  NOT  INCREHENT  KALI., 

DIMENSION  DAT(512)iBUF(256)»KBUF(256)  B  .2 
INTEGERI4  KRECilOSB 
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-►  DO  10  IDX*1>256 

f  KJUF  <  IDX)aOAT <  IDX>tB(JF<  IDX) 

|  DUF<IDX)*MT(IDX+2Vi) 

10  CONTINUE 
K*£C=KALL+1 

CALL  MtITEF(LUNfKR£CiKDUFf5l2>I0SB) 

RETURN 

END 
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SUBROUTINE  HILKIAT .PRES.RES.FRES) 

DIMENSION  DAT(256)fPR£S<15)»R£S(256)»FR£SU5) 

C - BY  BOB  DICK.  15  JAN  80.  UPDATED  16  JAN  80. 

C - RES  IS  RESULT.  ADJUST  ENDS  AS  FOLLONS. 

C - ABD  PRES  TO  END  15  POINTS  OF  PREVIOUS  RESULT. 

C - ADD  FRES  TO  FIRST  15  POINTS  ft-  NEXT  FUTURE  RESULT. 

C - 31  POINT  FIR  FILTER.  E1€N  TAPS  0.  8  DISTINCT  MAGNITUDES. 

C - GAIN  1HK-0.024  IN  0.03  TO  0.47  SAMPLING  RATE. 

DIMENSION  VARt316)»0UT<286).TAP<8) 

DATA  TAP/. 0205. .0213.. 0326.. 0488. 

2  .0730.. 1140.. 2040.. 6338/ 

->D0  10  IDX’1.30 
4  VAR(IBX)=0. 

|  VAR ( 28641 PX)=0. 

10  CONTINUE 

DO  20  IDX-1.256 
^  VAR(IDXt30)=DAT(IDX) 

20  CONTINUE 

DO  40  MID--16.301 
4  SUM-0. 

DO  30  IPT=1 .8 

4  SUN=SUN4( VAR(HID42IIPT-1 )- 
I  I  VAR<MID-2IIPT41))ITAP(9-IPT) 

30  CONTINUE 

0UT(MID-15)=SUN 
40  CONTINUE 
DO  50  IPX=1,15 
4  PRES(IDX)=OUT(IDX) 

|  FRrS<mX)=OUT(271+IDX) 

50  CONTINUE 

DO  60  IDX=1.256 
f  RES<  IDX  >  =OUT  <  15-FIDX  > 

60  CONTINUE 
RETURN 

END _ 

SUBROUTINE  LNPS(FM.DAT) 

DIMENSION  FM(6).DAT(512) 

C - LOUP  ASS  0.2  SAMPLING  RATE.  6  FOLE  PUTTERUORTH, 

C - FM  IS  FILTER  MEMORY.  SHOULD  BE  0  INITIALLY. 

DATA  AK1 . A1 1 , A21/ .27725. . 49575. - . 60474/ 

DATA  AK2 . A12 . A22/ . 20657 , , 36753 . - . 17582/ 

DATA  WC3.A13.A23/. 18007.. 32212. -.04240/ 

— ^DO  10  IPT=1»512 

4  YSAK1<DAT(IPT)4A11*FM(1)4A21!FM(2) 

U=Y42.*FM<mFM<2) 

FM(2)=FN(1) 

FM(1)=Y 

Y=AK2«  4 A1 2*FH ( 3 ) + A22IFH ( 4 ) 
y=Y42.*FMC3)4FM(4) 

FMt4)=FM(3) 

FM(3)=Y 

YaAK3IU+A13IFM(5)+A23IFM(6) 

y*Y42.IFN(5)4FM<6) 

FM(6)«FM(5) 

FM(5)*Y 
DAT(IPT)*N 
10  CONTINUE 

<— RETURN  B-14 
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- HTCQW  BY  BOB  DICK 

- FROM  NTSORT  28  AU6  80 

- COM  FILTERS  VOICE  DEPENDING  ON  NISTONING 

- DATA  Mtors: 

- CCR.CCIJ  LIN  NAG  SPECTRUM,  THEJN  COHPLEX  CORRELATION 

— wtin:  input  data  before  windowing 
— iwf:  input  buffer 

- KWF!  HEADER)  THEN  OUTPUT  BUFFER 

c — oruf:  data  overlap  OUTPUT  buffer 
C - TR.TI!  WINDOWED  DATA.  THEN  FFT.  THEN  FILTERED  DATA 

c — mm:  raised  cosine  <hannino  window 

BYTE  F1NAN(30),F2NAM<30) 

DINENSION  DAT  IN (512) .WNDW(5I2) >TR( 512) .11(5)2) 
DIMENSION  CCR(512),CCI(512),0BUF(254) 

DINENSION  JBUF(25A),KBUF(256),LRUF<25«) 

INTE6ERM  IREC, IOSB, IBF(9),KBFi?) 

EOUIVALENCE  (IBF.lBUF).(KBF.KBUf ) 

DATA  KBUF/254I0/DATIN/512I0./ 

DATA  0BUF/256I0 . /RAKE . CONB/ -0 .5.0.5/ 

DATA  PI, TWPI/3. 14159245, 6. 28318531/ 

DATA  SANPF9/6400./ 

C - GET  INPUT  FILE.  NUMBER'  OF  BLOCKS 

— >TYP£  10010 

ACCEPT  10020.LNTN.  (FlNAM(KH).KHsl  ,i.NTH) 
F1NAN(LNTNH)=0 

00P£N(UNIT=1.NAN£-F1NAN.TYPE='0LD'» 

1  BUFFERC01INT= - 1 .  ERR=9000 » READQNL  Y ; 

IREC=1 

CALL  READF ( 1 . IREC . I BUF .512. IOSB ) 

N1BLKS=IBF(9) 

TYPE  10030.N1BUS 
ACCEPT  10040, NHBLKS 

IF (NHBLKS.LT. 2) GO  TO  9000 - >ABC*T 

IF(NMBLKS.GT.N16LRS)G0  TO  9000 - >/ABO*T 

C - GF-T  OUTPUT  FILE,  WRITE  HEADER 

TYPE  10050 

ACCEPT  10020,  l.NTH,(F?NAN(KH),KH=l,LNTH) 

F2NAN(LNTHE1 1=0 

0CPEN(UNIT=2,NAME=F2NAM , TYPE* ' NEW ' , 

1  WJFFERCOUNT=-l,Efifi=?000) 

KBUF(1)=1 

KBUF(2)-6400 

SEC=NMBLKS/25. 

NIN=NH8LKS/1500 

ISEC=SEC*40.IMIN 

KIUF(8)-NIN 

KWJF(9)=ISEC 

K8F(9)=NrtBLKS 

CALL  WRITEF(2, IREC, KBUF.512, IOSB) 

C - SET  UP  DATA  WINDOW 

DO  10  IPX=1 ,512 

f  MNW(IDX)=0.5-0.5IC0S(PIIIDX/256.) 

10  CONTINUE 

C - GET  HISTUNIN6  TO  COMB  FOR 

TTPE  iOMO 
ACCEPT  10070, SHFFfl 

TYPE  10080  B-16 

ACCEPT  10070, PCTCH8 


FA GE  2 

CTAN8L«PI*PCTCMI/100. 

C - ZERO  OUTPUT  DUFFER 

M  15  IBXM.254 
f  KDUF<IDX>=0 
IS  CONTINUE 

C - READ  FIRST  DATA  RECORD  BEFORE  MAIN  LOOP, 

WIIN*1 

CALL  TFSADV(DATIN*NMIM) 

C - LOOP  MUMI£JLQF_RfCOR0$-l  TIMES, 

DO  80  NNIN=2,NNBLKS 

C - A--MAIM  LOOP.  DONE  ONCE  PER  DATA  RECORD. 

IKONON 

C - TAKE  TIME  HINDOOS  KITH  50  PERCE!*!  OVERLAP. 

CALL  TFSADV(DATIN»NMIN) 

DO  20  1DX=1»512 

A  TR(IDX)  'DATINdDX)tUNDW(iDX) 

I  TI(I0X)=O, 

20  CONTINUE 

C - FORM  AMPLITUDE  SPECTRUM 

CALL  CFFT<1*TR*TM> 

DO  30  IDX=1,?5*J 

A  CCR(IDX)=SORT(TR(IDX)*TR(IDXHTI(IDX)*TI(IDX)) 

I  CCI(IDX)=0. 

30  CONTINUE 

C - SET  NEGATIVE  FREOUENCT  PORTION  TO  ZERO 

DO  AO  IDX=257»512 
A  CCR(IDX)=0, 

I  CCI(IDX)=0. 

40  CONTINUE 

C - TAKE  INVERSE  FFT 

CALL  CFFT(-1»CCR»CCI*V) 

C - SEARCH  FOR  PEAK  FROM  LON  TO  HIGH  TIME  LIMIT 

PEAK=0. 

IDPK=25 

DO  50  IDX=25*80 

A  TENP<CR<I0X)*CCR<IDXHCCIUBX)*CCI(IDX> 

if<.not. <pem.lt. temphgq  to  45 

j  PEAK=TEMP 
f  IDPMIDX 
45  CONTINUE 
50  CONTINUE 

C - FIND  QUADRATIC  INTERPOLATION  PEAK, 

BS0=CCR(IDPK-1  )ICCR(IDPK-1  )+CCI  <  IDPK-  t)ICCI< IDPK-1 ) 
USIKCRI  IDPKEi  >ICCR<  IDPKE1 HCCI  ( JDPKtl  )*CCI  (IDPKE1 ) 
X*l. 

IF(P£AK,6T.US0)X*C. 

IF(BS0.GE,P£M)X*-1. 

IF(X,EQ,0,)X-0,5I(US0-BS0)/(2.*PEAK-BS<HJS0) 

TAU=IDPK-1+X 

C - INTERPOLATE  REAL  AND  I MAG  PARTS  (VIA  QUADRATIC) 

CFl=0,5t(CCR( IDPKE1 )-CCR(IDPK-l ) ) 
CF2=0,54(CCR(IDPK-1 ) ECCRl IDPKE1 ) ) -CCR( IDPK) 
PARTR*CCR<  IDPKHXNCF1EXKCF2) 

Cfl*0. 54(CCI  (IDPK+1 ) -CCKim-l ) ) 
CF2-0.5l<CCI(IDPK-l>ECCI(IDPKtl))-CCI(IDPK) 
jl  PARTI*CCKIDPK)«»<CFHX*CF2) 

PARTM*SORT  ( PARTRIPARTREPART I IPARTI )  B-17 
IF(PARTM,LE,0,)PARTM»1, 
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C - -COMPUTE  PHASE  ANGLES 

AN6LD*ATAN2(PARTI/PARTN,PttTR/PARTH> 
FRACA*TAU»SHFFQ/SAMPfG 
ANGLA^TGPII(FRACA-IFIX!FRACA) ) 
IF(ANGLA.GT.PI)ANGLA-ANGLA-fHPI 

C - MINUS  PI.LT.AW6LA.IE.PI 

DISTA*ANGLD-ANGLA 
IF  <  DISTA.LT .  0 .  )DISTA--fi  JSTA 
IF!DI$TA.6T.PI)DISTA-TGFTDISTA 
C - O.LE.DISTA.LE.PI 

C — . TEST  FOR  DATA  PHASE  GITH1N  CUT -ANGLE  Of  PREDICTED  PHASE 

SGITCH*COMB 

iF<  DISTA.CT  .CTAW6L  )SWITCH=RAKE 

C . CONI  OR  RAKE  ACCORDING  TC  PHASE  DIFFERENCE 

0NE6=PI*TAU/256. 

DO  70  IDXU.256 

A  CS=(PARTR*COS(OHEG<IDXHPARTIISIN(OHEG*IDX))/PARTM 
FILT=0.5+SUITCH$CS 
JDX'513-IDX 
TR( IDX)=FILT<TR( IDX ) 

TI<IBX.'-FILT*TT(IDX» 

TR( JDX)=FtLTtTR< JDX) 

Tl( JBX)=FJLT*TI ( JDX) 

70  CONTINUE 

C - TAKE  INVERSE  FFT 

CALL  CFFT(-l.TRiTI.9) 

C - OVERLAP  RESULTING  TINE  FUNCTIONS 

NM0UT=NNIN-1 

CALL  LAPOUT(TR.OBUF.?»NMQUT) 

80  CONTINUE 

C - WRITE  FINAL  DATA  RECORD  AFTER  MAIN  LOOP. 

DO  90  IDX= 1,256 

|  TR(IDX)=0. 

90  CONTINUE 

CALL  IAP0UWR.0BUF.2.NHBLKS) 

9000  CONTINUED - < - ABORT 

<—  STOP 

10010  FORMAT! '  HISTUNING  COMBER  INPUT  FJLF?') 

10020  FORMAT (Qi30A!) 

10030  FORMAT! lXi 18.'  BLOCKS.  HOG  MANY  DO  YOU  GANT?') 

10040  F0RMAKI8) 

10050  FORMAT ! '  GHAT  OUTPUT  FILE7' ) 

10060  FORMAT!'  GHAT  FREQUENCY  SHIFT?(USE  .)') 

10070  FORMAT <F10,2) 

10060  FORMAT!'  GHAT  PERCENT  COMBING?!USE  .)') 

EMD _ 

SUBROUTINE  TFSADV1X.ICALL ) 

C - BY  BOB  DICK.  REVISION  OF  14  MAR  80. 

C - ADVANCES  BUFFER  256  POINTS  PER  CALL. 

DIMENSION  X(512).IBUF<256) 

INTEGERS  IREC.IOSB 
->IREC*ICALLU 

CALL  READFU.IREC.IBUF.512.I0SB) 

DO  10  IDX* 1.256 
|  X(IBX)«X(IDX4?56) 

1  :«IDX4256)*IBUF(IDX) 

10  CONTINUE 
<— RETURN 
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SUBROUTINE  LAPOUT(DATrBUF>LUN>KALL) 

C - OUTPUTS  DATA  MITH  50  PERCENT  OVERLAP. 

C - INITIALLY  BUF  SHOULD  BE  ALL  ZERO. 

C - DOES  NOT  INCREMENT  UlL, 

OINENSION  DAT(512)fBUF(256>.KBlJF(250) 
INTE6ERM  KRECiIOSB 
— 10  IDX«1»25* 

T  KBUF<IDX)«DA1<IDX>fBUF<IDX) 

|  BUF< IDX) =DAT  ( I DX+256 ) 

10  CONTINUE 
KREC*KALL+1 

CALL  NRITEF(LUN.KREC»KBUF.512-I0SB) 

■♦“RETURN 

END 
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C - PtCES  BY  BOB  DICK 

C - FROM  PTHkS  5  SEP  60 

C - ESTIMATES  BOTH  PITCH  MD  HISTUNING 

C - INCREMENTS  HISTOGRAM  BY  CORRELATION  PEAK 

C - CAM  PLOT  2D  PITCH-NISTUNIN6  HISTOGRAM 

C--  ''ATA  ARRAYS! 

C - CCR.CCI*  LIN  MAG  SPECTRUM*  THEN  COMPLEX  CORRELATION 

r - DATINI  I'iPUT  DATA  BEFORE  HIMOONING 

C - HMT!  KiSTOGRAM  Of  MISTUNIN6  GUESSES 

C -  NFS  4tf0  TO  400  HZ  BY  5  HZ 

C - HNTPT!  HISTOGRAM  Of  HISTUNING.  PITCH 

C -  PITCH  18  BINS.  80  TO  260  HZ  BY  10  HZ 

C - RESERVE  2  PTS/IINE  FOR  PLOTTING  DATA 

C - .BUF1  INPUT  BUFFER 

C- . -TR.Ti;  NINDOHED  DATA.  THEN  FFT 

C— -  WNDU:  RAISED  COSINE  <1.  'RING)  HINDOO 
BYTE  F1NAH(30)<«NSH 

DIMENSION  DATIN  .  ir  i<512).TR(512).TI(512) 
DIMENSION  CCR(51?.»Ui(512)»HMT(14t ) 

DiMENfiON  IBUF(256).HMTPT(163»18) 

INTEGER**  IREC.I0SMBH9) 

EQUIVALENCE  (IBF.IBUf) 

DATA  DATI N/5 1 2*0. /T8PI/A . 2831853072/ 

DAT.T  HMT/ 1  u <0 .  /HNTPT/293410  >  / 

•-►TYPE  10010 

ACCEPT  10020.LNTH. (FINAM(KH) »KH*1 »LNTH) 
F1NAM(LNTHI1>=0 

00PEN<UNIT=1.NANE«F1NAM,TYPE*'GLD'» 

1  BUFFERCOUNT®- 1  *  ERR®9000 » READONLY ) 

IREC-1 

CALL  CEAUFU.tREC.IBUF.5l2.I0SB> 

N1BLKS=IBK(9) 

TYPE  i0030.NlPt.KS 


ACCEPT  10040. NFBLK 

JF(NFBLK.LE.0)00  TO  9000 - 

IFOtt-  KLK.GT .N1BLKS)G0  TO  9000 
ACCEPT  10040. NNBLKS 
NL  BL I , = NF  Bi ./  tNMBLK  S  - 1 
IF(NLBLK,Lr,NFBLK>GP  TO  9000‘ 
IF(NLBLK.G;.N1BLKS)GU  TO  9000 


> abort 

►ABORT 


■*.  ABORT 
►ABORT 


DO  10  IDX- i .512 

f  UNDH ( I r  X / ®0 . 5-0 • 5<C0S ( 3  *  1 4 1 59265* I DX/254 • ) 
10  CONTINUE 


SNPK®0. 

SMPT®0. 


S.*PT2=0. 

NFIRST^NFBLLU 


LAST*NLBLKM 
DO  80  NHIN=NFIiw. .LAST 
lREC*NKIN 

-TAKE  TINE  INDOHS  WITH  50  PERCENT  OVERLAP. 
CALL  TFSADV(DATIN.NNIN) 

20  IDXU.51? 

TR(IDX)eDATIN(IDX)*HNDH(IDX) 

1 t(IDX)*0. 

CONTINUE 

CALL  CFFTU.TR.7I.9)  13-20 

C-  — h— FORM  LINEAR  MAGNITUDE  SPECTRUM 


DO  80 

tlR 
-TA 

CA 

f 


20 
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f30  IDX-1  »2S6 

CCR(IBX)*SQRT(TR(IDX)*TR(1DX)FTI<IDX)*TI<IDX)) 
CCIUDX)«0. 

CONTINUE 

■SET  NEGATIVE  FREQUENCY  PORTION  TO  ZERO 
DO  40  !DX*257i5I2 
f  CCR(IDX)*0. 

CCI(IDX)*0. 

CONTINUE 

-TAKE  INVERSE  FFT 
CALL  CfFT<-l,CCR(CCI*9) 

PEAK«0, 

-SEARCH  FOR  PEAK  IN  8u  TO  254  T*RTZ 
DO  40  IDX*25>80 

,  TEHP*CCR<IDX)*CCR(IDX)+CCI(IDXXCCI(IDX) 

IF(. NOT, (PEAK, LT,1£HP))G0  TO  50 
j  PEAK* TEMP 
f  IDPK-1DX 
CONTINUE 
CONTINUE 

-FIND  QUADRATIC  INTERPOLATION  PEAK 
BSQ»CCR<IDPK-l)»CCR(IDPK-l)4CCl<IDPK-mCCI(IDPK-l) 
USQ-CCR( IDPKF1 )*CCR( I  DPMI HCCI < I  DPMI )*CCI ( !  DPMI > 
X»l. 

IF<PEAK,0T.IISQ)X*0. 

IF(ESQ.GE.PEAK)X*-1, 

IF(XiEQ.0t)X*0.5KUS2-88Q)/(2.tPEAK-BS8-USQ) 
TAU«IDPK-1+X 

-INTERPOLATE  REAL  AND  IHA6  PARTS 
CF 1 -0 . 5*  <  CCR  < IDTKF 1 ) -CCR ( I DPK- 1 ) ) 

CF2=0.5t<CCR<  IDPK-DFCCR!  IDPM1 )  >-CCR(  IDPK ) 
•PARTR*CCR(IDPK)M<CF1MCF2> 
CFl«0,5*<CCI(IDPK41)-CCmflPK-l)) 
CF2*0i5«CCI(ILPK-l)FCCI<IDPKFl))-CCI(IDPK) 
PARTI*CCI(IDPK)FXI(CFHXICF2) 

PARTH*SQR  T ( PARTRIPARTR  PPART I OPART I ) 

IF (PARTHiLEtOi )PARTH*1 . 

-FORK  INTERPOLATED  PEAK  FROM  PARTH. 

PEAK=PARTN*PARTN 
ShPK*SHPK4PEAK 

ANGL  =>AT  AN2  ( PART  I  /PARTH  > PARTR/PARTN ) 

PITCHc640Q./TAU 
SHPT*SHPT4PEAKIPITCh 
3HPT2=SHPT24PEAKIPITCHIPITCH 
SHF=PITCHIANGl/TMPI 
IXPT*(PITCH-70,)/10, 

1F(IXPT.LT.1)IXPT*1 
IF(1XPT.GM8)IXPT*18 
GO  TO  68 
CONTINUE 

J  SHT-SHF-PITCH 
CONTINUE 

IF(SHF.GE>-402.)GO  TO  65 
DO  70  IDX'l'li 
f  SHF*8HF+PITCH 

. ROUND  INSTEAD  OF  TRUNCATING 

V  IXHT-81.5P8HF/5.  b'Jl 

-IF (POINTER  OUT  OF  RAN6E)EXIT  LOOP 
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1F( .NOT. (IXMT.LE.141 ) )G0  TO  7S 
f  HAT < 1XMT )=HMT ( iXMT )+PTAK 
|  HrtTPT ( IXMT  * IXPT )=HMTPT (IXMT. IXPT) EPEAK 
70  CONTINUE 
75  CONTINUE 
80  CONTINUE 
AVPT-SHPT/SMPK 

SDPT =SQRT ( SHPT2/SHPK-AVPTIAVPT ) 

TYPE  1 0050  f  A1  ’°T .  SOPT 
PKMT=0. 
h  f  X = 8 1 

ini  100  IDX=1.141 

A  IF ( .NOT . (PKMT .1  T.HHTt inX) ) )G0  TO  90 
I  PXNTSHNT(IDX) 

\  «TX=IDX 
90  CONTINUE 
100  CONTINUE 

FRi1T=<«TX-81)$5, 

TYPE  10040. FKrtT 
TYPE  1007C 
ACCEPT  10080. ANS« 

IF(.NOT.(ANSM.EQ.'Y'))GO  TO  110 
I  CALL  PLTHST(HMTPT) 

^  ACCEPT  10080.ANSN 
110  CONTINUE 

9000  CONTINUED - * - AB6RT 

<-  STOP 

10010  FORMAT ( '  PITCH-HISTUNING  ESTIMATOR  INPUT  FILE?') 
10020  FORMATS.  30A1) 

10030  FORMAT < 1X» 18. '  BLOCKS.  GIVE  FIRST  <CR>  HOW  MANT.') 
10040  FORMAT  < 18) 

10050  FORMAT ( '  PITCH  AVG,  STD  DEV:'.2F8,1> 

10040  FORMAT! '  ESTIMATED  MISTUNIHG  (HZ)t'.F9.0) 

10070  FORMA'. i '  DO  YOU  WANT  A  GRAPH? <Y  FOR  YES)') 

10080  FORMAT (At) 


SUBROUTINE  TFSADV(X.ICALL) 

C . BY  BOB  DICK.  REVISION  OF  14  MAR  80. 

C - ADVANCES  BUFFER  256  POINTS  PER  CALL. 

DIMENSION  X(512).IBUF<254> 

INTEGERS  IREC.IOSB 
IREC=ICALlfl 

C ALL  RE ADF ( 1 . IREC . I BUT , 512 > IOSB ) 

DO  10  IDX=1 . 256 
A  ■  IDO=X(IDX+25A) 

|  X(IDXi?54)=IBUF(IDX) 

10  FONTINUE 
<"  RlTuRN 

END _ 

SUBROUTINE  PLTHST (HMTPT ) 

C - BY  BOB  DICK  4  AUG  80 

C - LAST  UPDATE  3  SEP  BO 

C - PLuTS  2D  HISTOGRAM  OF  MISTUNING  VS  PITCH 

DIMENSION  HMTPT (143. 18) » XCRD< 163) 

DIMENSION  AXISX(4).AXISY(4) 

DATA  AXISX/B1 . .81 , . 1 . . 13./AXISY/0. .9.2.0. >1 ./ 

C - FIND  HISTOGRAM  MAXIMUM  B-22 

->HMAX=0. 
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10 


DO  20  10X2*1  r  16 
DO  10  10X1-1 .161 
A  TEHP-HHTPT4 IDX1 > IDX2> 

|  IF(HHAX,LT,T£HP)HMX«TEHP 
CONTINUE 
20  CONTINUE 

—GENERATE  X  COORDINATES 
DO  30  1DX1«M61 
f  XCREKTOXn-IOXi 
3 0  CONTINUE 


XCRD<162)*1, 

C - 160  INTERVALS  AT  13/iN  IS  12,3  IN  HORIZ 

XCRD(163)*13. 

C - VERT  SCALE  IS  HIST  HAX  PER  0,8  INCH 

VSCALE-HHAX/0.8 

C - CLEAR  SCREEN)  RESET  ORIGIN 

CALL  PLOTS 


C~ 

C- 


C-- 


c 


DO  40  LNE®1)18 
-if-DONE  ONCE  PER  LINE 

. 17  INTERVALS  AT  0,5  INCH  EACH  IS  B.5  IN  VERT 

TORO*(18-LNE)tO,5 
CALL  PLQT(0,iY0RG»-3) 

HHTPT(162.LNE)*0, 

HHTF'T  ( 163»LNE )  sVSCAL£ 

. PLOT  LINE 

CALI  UNE(XLRD)HHTPT(1)LNE>)  161fl»0)0> 
YORG—YORG 

CALL  PLOT (0, ) VORG) -3) 

40  CONTINUE 

-PLOT  0  HZ  HISTUNING  AXIS 


CALL  LINE(AXISX)AXISY)2) 1)0)0) 

C - REPEAT  LAST  PLOT  LINE  TO  FILL  BUFFFR 

CALL  LIN£<XCRD)HHTPT<1)18M61)  1)0)0) 
<-  RETURN 
END 


B-23 
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C - VCRfSP  BY  BOB  DICK  2?  NAY  80,  LAS!  UPDATE  30  HAY  80. 

C - FROH  VCRFPL  29  HAY  80, 

C - VOICE  RE -FILTER  AND  SEPARATE. 

C - LATTICE  FILTER,  THRESHOLD  RESIDUAL,  RECONSTRUCT 

C - DATA  ARRAYSi 

c — brsd:  backward  residual 

C - COEFi  LATTICE  FIl TER  COEFFICIENTS 

c — frsd:  forward  residual 

c - RFCOEF!  RECONSTRUCTION  FILTER  COEFFICIENTS 

c — sigdif:  input  signal  afiek  differencing 
c — sigin;  input  signal 

C - SIGOUT!  OUTPUT  (RECONSTRUCTED)  signal 

C - TRSD!  FORWARD  RESIDUAL  AFTER  THRESHOLDING 

DIHENS  ION  J»RSD<  1 50 ) ,  COEF  ( 20 ) ,  FRSD  U  50) ,  RFCOEF  i  it ) 
DIMENSION  SIGDIFI 150) ,S1GIN<  150)  ,SIGOUT(  150) , TT(SD(  150) 
DIMENSION  BLANKS ( 64 ), SI GSEC< 150) 

INTEGERS  I  RE C ,  I  OSS ,  I BF  <  9 ) ,  KBF  ( 9 ) 

BYTE  F1NAH(30),F2NAH(30),F3NAH(30) 

DIMENSION  18UFI256) »KBUF(25A) 

EQUIVALENCE  (JBF,IBUF),(KBF,KBUF) 

DATA  KBUF/256I0/BLANKS/64I0./ 

— ►TYPE  10010 

ACCEPT  10020 » LNTH ,<  r I NAH ( KH ), KH' 1 , LNTH , 

FINAHILNTH+D^O 

OOPENUWIM ,NAM£=F1NAM»  TYPE= ' OLD ' , 

1  BUFFERC0UNT=-1 , ERR- VOOO, READONLY/ 

1REC*J 

CALL  READF(|,!REC,iDU),5i:,I0SB) 

N1»LKS'IBF<9) 

TYPE  10030, N1BI.KS 

ACCEPT  10040, NHDLKS 

IF<NHBLKS.LE.0)NN8LKSsNlBLKS 

IF<NHBLKS,0E,N1DLKS)NHBLKS*N1BLKS-I 

NHLPS*NHBLKSI? 

TYPE  10050 

ACCEPT  10020, LNTH,  <F2NAH(KH),KH=1  ,l.NTH) 

F2NAHUNTHI1  )=0 

OOPEN  ( UNI  T«2 ,  NAHESF  2NAH ,  T  YPF.= '  NEW ' , 

1  BUFFERC0UNT--1 ,ERRS9000) 

ACCEPT  10020, LNTH,  (F3NAH(KH),KHM, LNTH) 

F3NAHIINTHM  )=0 

00PEN(UNIT-8,NAH£ZF3NAN, TVPF-'NEW' , 

1  BUFFERCOUNT*-l,ERR-VtOOi 

KDUF(1)-I 
KBUF(2)-6400 
SEC*NHBLKS/25, 

MIN-NHBLKS/IGOO 
ISFC=SEC-60.*MIN 
KBUP(8)-HIN 
KBUF ( 9 ) = i SEC 
KBFI9)-NHBLKS 

CM.L  WfITEF(2,IREC, KBUF, 512,1088) 

CALL  WPITEFIP-IREC, KBUF, 512, IOSB) 

TY,E  10060 
ACCEPT  10070, RSTH 
CALI  INTGD(148,?0) 

CALL  WROUT (BLANKS, 20,1,2)  B-24 

CALL  WR0UT(8LANKS,20,7,8) 


20 


30 


30 


-THIS  LINES  UP  OUTPUTS  WITH  INPUTS. 

PO  40  KLPS-l.NNLPS 

CALI.  8£TMT<  1  f SIOIMf ICOCC) 

-I CODE  NEGATIVE  MEANS  HO  HOPE  DATA 
IF!.NOT.!ICODE.GE.0))O0  TO  70 

C - ^--DIFFERENCE  SIGNAL  FOR  PRE-EHPHASIS 

SI6DIF<1)«SI6INU> 

DO  7.0  IDX«2rl48 

f  S(GDIF(IDX)»SIGIN(IDX)-SIGIN(IPX-1) 

CONTINUE 

CALL  LTFl.T<SlGDJF»14B»COfF»20fFRSD»BRSB) 

CALL  8LNKR!FRSD»TRSD.  t48.RSTH) 

CALL  RL  TFL  T ( TRSP . 1 48 . COEF . 20 . S IGOUT ) 

C - +--SUH  SIGNAL  TO  UNDO  PRE-EHPHASIS 

DO  30  JDX-2f 148 

f  SI GOUT ( IDX)*SI60UT ( IDX) I SIGOUT ( IDX-1 ) 

CONTINUE 

CALL  NR0UT(SJG0UT(21)»128il»2) 

DO  30  IDX«21»H8 

f  SIGSEC(I0X)*SIGIN(IDX)-SIG0UT(1DX) 

CONTINUE 

CALL  WROUT <SI6SEC< 21 )»128>2»8) 

60  CONTINUE 
70  CONTINUE 
9000  CONTINUE 
-NE — STOP* 

10010  FORMAT!'  INPUT  FILE  FOR  VOICE  RECONS TRUCTQR7' ) 

10020  F0RMAHQ.30A1 ) 

10030  FORMAT < IX .18 » '  BLOCKS.  HON  MANY  DO  YOU  HANTT') 

10040  FORMAT! 18) 

10050  FORMAT!'  HHAT  TWO  OUTPUT  FILES?') 

10060  FORMAT!'  MULT  OF  STD  Df.V  FOR  RESID  THRESH?! USE  .)') 
10070  FORMAT(F12.5) 

_ END _ _ _ 


C~ 

C- 

C~ 

C- 

C- 

C— 

c- 

c— 

c- 


SUBROUTIHE  GETDAT ( LUN . RDAT . I CODE ) 

-BY  DOB  DICK  1  MAY  80, 

-READS  OVERLAPPING  DATA  BLOCKS, 

-IBUF  IS  INPUT  BUFFER.  POINTER  KIB 
-DBUF  IS  DATA  BUFFER.  POINTER  KDB 
-RDAT  IS  DATA  RETURNED 
-NSIZ.MISZ  IS  DATA  SIZE  PER  CALL 
-NVLP.MVLP  IS  DATA  OVERLAP  BETWEEN  CALLS 
-INTGD  IS  ENTRY  TO  INITIALIZE  AND  SET 
- SIZE  AND  OVERLAP. 

DIMENSION  IBUF! 256) .DBUF! 1024) .RDAT! 1024). I Q5B12) 
BYTE  ERRBYT 

EQUIVALENCE  (ERRBYT. IOSB) 

INTEGERI4  KREC 
>NFIRST=KDB+1 
DO  30  KDB«NFIRST.MSIZ 

IF(,NOT,!KIB,EQ.O))GO  TO  10 
4— HERE  NEED  MORE  DATA 
KR£C*KR£CH 

CALL  READEa.UN>KREC»I&lf>512.I0SB) 

KIB-1 

IF!ERRBTT.Ut.,0)KIB»l 

CONTINUE  B-25 

IF(,NOT.!KIB.GT.O))60  TO  20 
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C . +--HERE  THERE  IS  NO  END  Of  FILE 

DBUF<KDB)*IBUF(KIB) 

I  KIB-KIBf 1 
Y  JF<KIB.6T.254)KIR*0 
20  CONTINUE 
30  CONTINUE 

IF(. NOT.OUB.8E. OHIO  TO  /O 
C— HERE  NO  EOF 

DO  40  KRD*1»HSJZ 
|  RDAT(KRD)=DBUF(KRD) 

40  CONTINUE 

IF(.NOT.(HVLP.GT.O))GO  10  40 

C . -j— HERE  THERE  IS  OVERLAP 

DO  50  KDB=liHW.P 

|  DBUF I KDB ) =DBUf ( KDBFHS1 Z  -HVl.P ) 
50  y  CONTINUE 
40  CONTINUE 
'  '  KOB=MVLP 
70  CONTINUE 
ICODEaKlB 

IF(ICQDE.GT.O)ICODF-C 
<—  RETURN 

ENTRY  INTGD(NSIZiNVLP) 

C - INITIAL  ENTRY 

->HSIZ*1 

IF<1.LE.NSIZ.AND.NSJZ,IE.1024)HSI7*NSIZ 

HVLP=0 

IF  ( 0 .  LE .  NVLP .  AND .  NVI.  P .  l.T .  NS  1 1 )  HVLP*NVIP 

C . SKIP  HEADER  BLOCK 

KREC-1 

KIB»0 

K0B*0 

RETURN 


SUBROUTINE  LTFLT(SIGiLENiCOEF)NORDER)F!?SDfPRSD 
-BY  BOB  DICK  1  HAY  80.  LAST  UPDATE  2  HAY  80 
--NEH  HORE  GENERAL  VERSION  OF  ADFl.T . 

--SIU  IS  INPUT  UAVEFORH  LENGTH  LEN 
--COEF  NHL  BE  LATTICE  FII.TER  COEFS.  NORDFR  OF  THEH 
--FRSO  HILL  BF  FORNARD  RESIDUAL 
--BRSD  MILL  BE  BACKWARD  RESIDUAL  (FIRST  NORDER  ZEROED) 
-DOES  NOT  FTSELF  DIFFERENCE  SIGNAL  '  OR  PRE-EMPHASIS 
DIHENSION  SIG( 1024)iC0EF(32) iFRSIU  1024) (BRSDl 1024 ) 

►  IF ( .NOT • ( 1 . LE . NORDER . AND . NORDER . LE , 32  /)  GO  10  9000 — 
IF< .MOT • (HORDER+1 .LE.IEN.AND.IEN.LE.1024))G0  TO  9000- 
FRSD(1)=SIG(1) 

BRSD<1)=0. 

flO  IDX=2iLEN 

FRSD(IDX)*SJG(IDX) 

BRSD<IDX)*SIG(I0X-1) 

0  CONTINUE 
»  40  I TER=1 .NORDER 
SFS0*0. 

SBSfcO. 

DOTPR*0, 

DO  20  IDXMTERH.LEN 

^  SFSQ=SFSOEFRSD(IDX)»FRSD(IDX)  B-26 
SBSO^SDSOEBRSD ( IDX  > I BRSD ( IDX ) 


A&OAT 
-:J►  ABORT 


K**  ■ 


VCRFSP.FTNU 


t  DOTPR*DOTPR4FRSD<IDX)*BRSOUDX) 

CONTINUE 

nENON-SFSAfSBSO 

IF(DENOH,LE.O.)DENOH»1. 

COEF ( I IER ) *-2 . IBOTPR/DENOH 
WOLB^O, 

DO  3-)  IDX=lTEfi+lfLEN 
A  BRTHP*COEF (ITER) IFRSB ( I DX ) +  RRSD 1 1 DX ) 
FRSD( IDX)-FRSB(  IDX ) FCOEF ( 1 TER)IBRSD( IDX) 
BRSD(IDX)«BROLD 
BROLD'BRTHP 
CONTINUE 


40  CONTINUE 
9000  CONTINUE- 
-4— RETURN 


- ABORT 


SUBROUTINE  BLNKR(DATINtDATQUT.LENiTHRES) 

C - BT  BOB  OICK  1  NAT  80 

C- . DATOUT»DATIN  WITH  SHALL  VALUES  ZEROED. 

C - VALUE  IS  SHALL  IF  MITHIN  STD  DEVITHRES  OF  ZERO 

C - DAT  IN  CAN  ALSO  BE  WATOUT. 

DIMENSION  DATINt 1024) i DATOUT l 1024) 

~>IF(.N0T.U.LE.LEN.AND,LEN.LE.10:4))00  TO  9000 - $►  ABORT 

SUHSQ*0. 

DO  10  IDX-lrLEN 

|  SUHSO=SUHSO+DATINaDX)IDATTH(ITiX) 

10  CONTINUE 

STDV-SORTISUHSO/LEN) 

CUTF*THRESBSTDV 
DO  20  1BX*1.LEN 
A  RESULTS. 

IF(DATIN(1DX).LE.-CUTF)RESULT*DATIN(IDX) 

IFICUTF.LE  .DATIN(  IDX))RESUI.TsDATINi  IDX) 

DATOUT< JDX)*RESULI 
20  CONTINUE 

9000  CONTINUE "4 - ABORT 

RETURN 

_ END _ 

SUBROUTINE  RLTFLT(DATiL£N»CQEF»NORDER»RES) 

C - BY  BOB  DICK  2  HAY  80 

C - DOES  RECIPROCAL  OF  LATTICE  FILTERING 

C - DAT  IS  INPUT  (APPROX  OF  LATTICE  FORWARD  RESIDUAL) 

C - LEN  IS  LENGTH  OF  INPUT  AND  OUTPUT 

C - COEF  IS  VECTOR  OF  LATTICE  FILTER  COEFFICIENTS 

C - NORDEK  IS  NUHBER  OF  FILTER  COEFFICIENTS 

C - RES  IS  OUTPUT 

C - RRSD  IS  REPRODUCTION  OF  VECTOR  OF  BACKWARD  RESIDUALS 

C - DOES  NOT  SUN  OUTPUT  TO  UNDO  PRE -EMPHASIS 

DIMENSION  DAT(J024bC0EF<32)»RES<1024>,BRSD<33> 

— >IF(.NOT.<t.LE.NORDER.AND.NORDER,LE.32))GO  TO  9000 - ►ABORT 

IF( .N0T.IN0RDERF1.LE.LEN.AND.LEN.LE.  1024))GO  TO  9000 - ►ABORT 

X!  30  KDATM.LEN 
'k  FNEX*DAT<KDAT) 

!F<.M0MKB*T,GE.2))6C  TO  20 
NSTG-KDAT-I 

IF(NSTE.GT.NORD£R)NSTG*NORDER 
DO  10  JST6sl.»NSTG  B-27 

h  A  KST6«NST&H-JSTG 


A 
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10 

20 


A 

FPR£*fN£X-COEF<KSTG)<BRSO<KSTG) 

BRSD<KSTG+l)=COEF(KSTfi)*FPRERBRSD(KSTG) 

FN£X=FPR£ 

’’  CONTINUE 
CONTINUE 
RES<KMT)ZFN£X 


I  BRSD(1)=FNEX 
30  CONTINUE 

9000  CONTINUED - < - ABORT 

-RETURN 

ENP _ 

SUBROUTINE  WROUT(PATfNUM>ICHANfLUN) 

C - BY  BOB  DICK  29  HAt  80.  LAST  UPMTE  30  NAT  80. 

C - MULTICHANNEL  OUTPUT  BUFFERER  AMD  WRITER 

C - DAT  IS  OUTPUT  DATA 

C - NUN  IS  NUMBER  OF  POINTS 

C - tCHAN  TELLS  WHICH  OF  4  CHANNELS 

C - LUN  IS  OUTPUT  LOGICAL  UNIT  NUMBER 

C - INWRT  IS  REINITIALIZATION  ENTRY 

INTEGERS'!  JOSB.KfiECM) 


DIMENSION  DAT<768MRUFU024iA»KBUF(4> 


DATA  KBUP»KREC/4*0.4*2/ 

“>IF(.N0T.(1.LE.NUM.AN0.NUM.LE.1024)/G0  TO  9000 - »-ABO*.Tl 

IF(.N0T.(1.LE.ICHAN.AND.ICHAN.LE.4>)G0  TO  9000 - >ABORT  t 

LBUf*KBUF(ICHAN) 

DO  10  KDAT-liNUM 
^  L8UF=L8UPT1 
I  IBUF(LBUFfICHAN)=DAT(KDAT) 

10  CONTINUE 
20  CONTINUE 

IF(.N0T.<LBUF.GE.258))G0  TO  40 

{CALL  WRI TEF (LUN>KREC( ICHAN) . IBUF( 1 i ICHAN) * 512* IOSB) 
KR£C(ICHAN)=KREC(ICHAN)E1 
DO  30  IDX=l»76fl 


30 


f  IBtif  (IDXr  ICHAN)-JBUF(IDXT256rICHAN) 
CONTINUE 


y  LBUFHBUF-25A 
GO  TO  20 
40  CONTINUE 


KBUF(ICHAN)=t.6Uf 

9000  CONTINUE  - ABORT  1 

RETURN 

ENTRY  INWRT(ICHAN) 

->IF<.N0T.a,LE.ICHAN,AND.ICHAN.lE.4))G0  TO  9100 - >AB0RT2 

KBUF(ICHAN)=0 


KREC(ICHAN>=2 

9100  CONTINUE^ - ABORT Z 

RETURN 

END 
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C - NISTUNED  SSI  MAIN.  LAST  UPDATE  6  OCT  80 

C - WITTEN  21  MAR  80 

C - DATA  ARRAYS! 

C - DBtR.DBLI!  DOUSLEl)  SAMPLE  RATE  REAL.  IHA6  PARTS 

c — fre.fih:  future  data  real,  ihag  parts 

C - FMR.FMI5  LOWASS  FILTER  MEMORIES  FOR  REAL.  IHAG  DATA 

C - FMO!  LOWASS  FILTER  MEMORY  FOR  OUTPUT  DATA 

C - FRES!  HILBERT  FILTER  FUTURE  RESULT  (TRAILER) 

c — iuf:  input-output  buffer 
c — olre.olih:  old  DATA  REM.)  ihag  parts 

C - PRES!  HILBERT  FILTER  PREVIOUS  RESULT  (LEADER) 

C - REO:  (REAL)  OUTPUT  DATA 

BYTE  6MAME(30).FNAHE(30) 

INTEGERS  IREC.I0SB.IB(9) 

DIMENSION  I8UF(256).FR£(25A)»FIM<254> 

DIMENSION  PRES(15)»R£S(25A)>FR£S(15) 

DIMENSION  0LRE(256)»0LIN(254).RE0<256) 

DIMENSION  rMR<6).FMI<6).FH0(6)iDBLR(S12)iDK.I(SI2) 
EQUIVALENCE  (IB.IBUF) 

DATA  FR£.FIHiFR£S/254t0..256tO>.15tO>/ 

DATA  FMR.FMI .FMO/6YO . .610. »AIO./ 

C - GET  DATA  FILE  NAME 

->TYPE  10000 

ACCEPT  10010.LNTH.(FNAME(KH).RH=1.LNTH) 

FNAME(LNTH+1)=0 

C - OPEN  UAVEFORM  FILE.  READ  SIZE. 

OPEN( UNI T=2 ,NAM£=FNAME . TYP£= 'OLD ' . 

2  BUFFERCOUNT = - 1 . ERR=9000 .READONLY ) 

IREC=1 

CALL  READf (?>1R£C.IBUFi512iIQSB) 

NMBLKS=IB(9) 

TYPE  10013.NMBLKS 

C - GET  NUMBER  OF  BLOCKS  TO  PROCESS. 

ACCEPT  10017.NBLK 

IF(.N0T.(NBLK.GE.1))G0  TO  9000 - 9».A80*.T 

IF(.NOT,(NBLK.L£.NMBLKS))GO  TO  9000 - >-ABOR.T 

C - GET  MISTUNING  TO  SIMULATE. 

TYPE  10020 

ACCEPT  10030.FQSF 

0MEGA=3. 141592G5AIFOSF/6AOO. 

AN6LE=0. 

C - OPEN  OUTPUT  FILE.  WKITF  HEADER. 

TYPE  10090 

ACCEPT  JOOIO.LNTH. (6NAME(KH) »KH=1 .LNTH) 

GNANE(LNTHil)=0 

OPEN(UNITsl.NAME=GNAME.TYPE='NEU' . 

2  BUFFER C0UNT=-1.ERRS9000) 

DO  10  IDX=1.256 
f  I8UF(IDX)=0 
10  CONTINUE 
IBUF(1)=1 
IBUF(2)=6400 
SEC=NBLK/2S. 

HDMOLK/1500 

IBUF(8)=MIN 

IIUF(9)=SEC-60.tMIN 

IB(9)*NBLK  B-29 

CALL  WITEFU.IREC.IBUF.512.I0SB) 
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T*I*2.*3, 141592654 

C - LOOP  NUMBER. OF.  H.0CXS41  TINES. 

LAST=NBLK42 
DO  ISO  JBLK=2,LAST 

C - ^ — IF(WOT  LAST  TIHE)6£T  BATA.  ELSE  USE  ZEROES, 

IF ( .MOT . < JBLK . LT .LAST ) >60  TO  JO 
-HERE  HORE  DATA.  USE  IT. 


30 


40 

50 


C— 


60 


C— 


70 


90 


C~ 

C— 


100 


IREC=J8U 

CALL  R£ADf (2,IR£C,IBUE,512,10S8> 

SO  TO  SO 
CONTINUE 

-HERE  NO  NQR£  DATA.  USE  ZEROES. 

DO  40  IDX--1.256 
IfcUf (IDX)=0 
CONTINUE 
CONTINUE 

— AGE  THE  DATA  (OLD: FUTURE i  FUTURE .*=«£«) 

DO  60  I OX  s I r  256 
OLRt(IDX)=FR£( IDX) 

0LIH(IDX)-FI«(IDX) 

FR£(IDX)=IBUF(1DX/ 

FI«(IDX)-0. 

CONTINUE 

--SAVE  PREVIOUS  HILBERT  FILTERING  TRAILER. 

BO  70  IDX=1,15 
f  FIH(1BX)=-FR£S(IDX) 

CONTINUE 

--DO  HILBERT  FILTERING. 

CALL  HILB(FRE>PRES,RES,FRES) 

DO  80  IDX= 1,256 
f  FIM(IDX)-FIN(IDX)-RES(IBX) 

CONTINUE 

-UPDATE  PREDICTION  (HI. BERT  FILTERING  LEADER). 
BO  90  IDX-1 i 15 

f  OLIN( IOX+241 )=OLIM( IDXP241 ) -PR£S( IDX ) 
CONTINUE 

-DO  TWO  INPUTS  BEFORE  FIRST  OUTPUT. 

-IF(FIRST  TINE  IHRUJSMP  OUTPUT. 
IF(.HOT.UDLK.GT.?))GO  TO  140 

■HERE  NOT  FIRST  BLOCK,  PRODUCE  OUTPUT. 


-DOUBLE  SAMPLING  RATE  BY  INTERPOLATING  ZEROES. 
DO  100  IDX-1. 256 

DBLR( IDXI2- 1 ) sOLRE (IDX; 

DBLR(IDXl2)-0. 

DBLKIDXI2-I  )=OLIrt(  IDX) 

DBLI(IDXI2)r0, 

CONTINUE 

SMOOTH  VIA  LGNPASS  FILTER'. 


CALL  L«PS(FNR,DBLR) 

CAU  LyPS(FNIiDBLI) 

C - -|— j — MULTIPLY  BY  COMPLEX  EXPONENTIAL,  KEEP  REAL  PART. 

DO  110  1DX=1,512 
AN6L£=AN6LEP0M£GA 
IF( AN6LE ,GT . TNPI ) ANGLE-AN6LE- TUPI 
IE(ANGLE.LT.-iyPI)AN6L£=ANGlEf TUPI 
D BLR( IDX) =DBLR( IDX) IC05( ANGLE i- 
OBLKIDXUSINt  ANGLE) 

110 |  |  CONTINUE 

B-30 
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< 
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C - 


C_.j_.x_, 


--LONPASS  TO  PREVENT  ALIASING. 


CALL  LUPS(FHOiONJt) 

— | — HALVE  THE  SAMPLING  RATE. 

DO  120  19X-1 >?S6 
f  REO< I OX ) =DBlR< IDX82 ) 

CONTINUE 
ITE  RESULT. 

DO  130  MX=1»?54 

-NULT  FT  ?  TO  COMPENSATE  FOR  ATTEN 
IWJF(I0X)=2.IRE0(IDX) 

130 1  |  CONTINUE 

IR£C=JBU-i 


HJ 

t 


CALL  URITEFdilRECiIDUF .512.  JOSfD 
140 1  ’CONTINUE 
150  CONTINUE 
CL0SE!UNIT=2) 

CLOSEOJNITd) 

TYFE  10050 

9000  CONTINUE < - ftBOP.T 

<—  STOP 

10000  FORMAT ( '  FILE  NAME  FOR  SSB  SIMULATION?') 
10010  F0RMAT(Q»30A1) 


10013  FORMAT! '  HAVE', 18.'  ROCKS.  HON  MANY  00  YOU  «ANT?') 
10017  F0RMAT(I8) 

10020  FORMAT!'  FREQUENCY  SHIFT  OF  HISTUNING?(HZ)') 

10030  F0RMAT(F8.0) 

10040  FORMAT! '  FILE  NAME  FOR  OUTPUT?') 

10050  FORMAT!'  FILE  HRITTEN. ') 

END 
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C - WTMEKG  BY  BOB  DICK  11  MR  80.  UPDATED  IS  APR  80. 

C - DOES  WEIGHTED  MERGE  OF  TWO  FILES. 

BYTE  F 1 HAH <  30 ) »F2NAM(30) ,F3NAM<30) 

DIMENSION  IBUF(256)»JBUF(236),KBUF(256) 

INTEGERI4  IREC»I0SB»IDBF<9),JDBF(9>,KDIF<9) 
EQUIVALENCE  (IDBF,IBUFMJDBF,JBUF),(KDBF,KBUF) 
DATA  F I NAN/30I0/F 2NAM/.10I0/F3NAM/30I0/ 

DATA  KBUF/25AI0/ 

— >-  TYPE  10010 

ACCEPT  10020,F1NAN 
F1NAN(30)=0 

00F'EN(UNIT=1,NANE=F1NAM»  TYPE* 'OLD' » 

1  BUFFERC0UNT*-l»£RR*9OOO. READONLY) 

ACCEPT  10020. F2NAM 
F2NAIK30U0 

OOPEN ; UNI T=2 , NANE=F2NAM, TYPE* ' OLD ' » 

1  BUFFERCQUNT*- 1 . £RR*9000 » READONL  Y ) 


IREC--1 

CALL  READF(1,IREC,IBUF .512.I0SB) 
CALI.  READF (2.  IREC.JBUF.512r  IOSB) 
N1PLKS=IDBF(9) 

N2BLKS- JPBF  <9 ) 

NSPLKS-N1PLKS 

IF(NSELKS.GT.N2BLKS)NSBLKS*N2BLKS 
TYPE  10030, NSBLKS 
ACCEPT  10040, NOBLKS 
I F ( NOBLKS . LE . 0 ) N0BLKS=NS8LKS 
I F ( NOBLKS . GT . NSBLKS ) NOBLKS*NSBLKS 
RMS 1=0. 

RMS2--0. 


DO  20  IBLK=1, NOBLKS 
A  IRECSIBLKF1 

CALL  READF(1,IR£C,IBUF»512,I0SB) 
CALL  READF(2»IREC. JBUF.512.I0SB) 
BSUM1=0. 

BSUM2=0. 


DO  10  IDX=1,256 
X=IBUF(IDX) 
BSUN1*BSUM1FXIX 
X=JBUF(IDX> 
BSUN2=BSUM2FXIX 
10  |  CONTINUE 

RHS1=RHS1 +8SUN1/256 . 
RHS2=RMS24BSUM2/256. 

20  CONTINUE 

RMS1=SORT(RMS1/NOBLKS) 

RH52=SQR  T  <  RMS2/N0BLKS ) 

PB2*20  ,ILOC(  RHS2/RMS1 )  /LOG  ( 10 , ) 
TYFE  10050.RMS1.RMS2.DB2 


TYPE  10060 
ACCEPT  10070, GNDB 
TYPE  10080 
ACCEPT  10020.F3NAM 
F3NAM<30)*0 

IF(.N0T.(F3NAM(1).NE.'  '))60  TO  9000 - >-/RBORT 

0CPEN(UNIT*B,NAME*F3NAM» TYP€*'NEU' » 

1  BUFFERCOUNT*-1»ERR»9000)  B-32 


KBUF(1)=1 


MUF<2)»4400 

SEMI0BIKS/23. 

HUMNOBlKS/1500 

ISEC-SEC-AO.MIN 

KNIF\8)*MIN 

KBI)F<9)=IS£C 

KDBF<9>-H0BLKS 

IREC=1 

CALL  WRITEF(0rIR£C»KBUF>512»IOSB) 
PAIN»10.M(GNDB/20.) 

FCTR1=1 ./(l.TGAIN) 

FCTR2»GAINtFCTRl 
DO  40  I8LK*1»N0BLKS 
IREC=IBLK+1 

CALL  READF(lfIREC»IBUF»M2fIDSB) 

CALL  READF  <  2 » IREC » JBUF  >  512  > 1088 ) 

DO  30  IDX-li254 

4  K8MF(  IDX )  *FCTR1IIBUF(  IDX )  TFCTR2IJBUF  <  IDX ) 

30  CONTINUE 

CALL  NRITEF(8>1REC»KBUF»512»IQSB) 

40  CONTINUE 

TYPE  10090>NOBLKS 

9000  CONTINUE < — - ABORT 

4—  STOP 

10010  FORMAT  < '  TWO  FILES  FOR  WEIGHTED  MERGE? ' ) 

10020  F0RNATT30A1) 

10030  FORMAT < '  NAME ' »I8» '  BLOCKS.  HON  HANY  TO  MERGE? ' ) 
10040  FORMAT (18) 

10050  FORMAT ( '  RMS'»2F10.2. 2N0/1ST1 '.F10.2> '  DB.') 
10040  FORMAT ( '  MHAT  DB  GAIN  FOR  2ND  FILE?') 

10070  F0RHAT(F10,2> 

10080  FORMAT < y  WHAT  OUTPUT  FILE?') 

10090  FORMAT < '  HEADER  AND'iI8>'  DATA  BLOCKS  HRITTEN.') 
END 
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IC - VTNHR6  ir  Ml  IIOC 

C - FROM  OTHERS  29  MB  M 

C - MES  ICI6KTO  ICR6E  OF  TM  FILES  MB  MIS  GAUSSIAN  NOISE. 

OTE  F1NAM(30),F2NAN(30),F3NAH(30) 

l  ‘  DMCNSJQN  IMF(236),JMF(2S6),KMT<256) 

I  INTE6ERI4  IR£C,I0$I,IMF<9),JUF(9),KISF<9> 

f  EMJ1  VALENCE  <IMF,IIUF),(JllF,JWFM»»f ,KJUF) 

t  DATA  F lKAH/30tO/F2NAH/30tO/F3NAN/30IO/ 

|  DATA  K8UF/256I0/ 

[  ->  type  10010 

ACCEPT  10020,  FINAfl 
F1NAH(30)=0 

OOP£N(UNIT=l»NArtE=FlNAH,TYPE='OLD'> 

1  BUFFERC0UNTs~l, ERR*9000, READONLY) 

ACCEPT  10020, F2NAH 
F2NAH(30)=0 

00PEN(UNIT=2,NA*=F2MH,TYPE='0LD' , 

1  8tFf£RC0UWT=-l,ERR=9000,READ0dr) 

IREC=1 

CAU  R£ADF(1»IR£C»IBUF,512»I0SI) 

CALI  READF(2»IREC,JIUF»512»I0SD) 

NlMJCS=IDBf  (V) 

N2DLKS=JD8F<9) 

NSBLKS-NIRLKS 

IF  (NSBLKS  ,6T  .N2BLKS)NSKJCS=N2BUCS 
TYPE  10030, NSBLKS 
ACCEPT  10040, N0BLKS 

j  IF(N0BLKS.LE.0)N0BLKS=NSBLKS 

I F ( M0BLKS . GT . NSBLKS ) N08LKS=NSBLKS 
RHS1=0. 

RHS2=0. 

i.  DO  20  IBLK=1,N0RLKS 

IREC=IBLKE1 

CALL  READR1,IR£C,IBUF,512,I0SI) 

CALL  fc£ADF<2,IREC,JKJF,512,I0SI) 

BSUH1=0. 

BSUN2=C. 

DO  10  IHX=1,256 
^  X=IBUF(IBX) 

DSUHlsBSWU+XIX 
X=JBUF(IDX) 

BSUM2=BSUH2+XtX 
CONTINUE 

RHSJ'RNSl+BSUHl/256i 
RftS2=RAS24BSUH2/256, 

20  CONTINUE 

RH$l=SfiRT(RH$l/NOILKS) 

RHS2=SQRT(RHS2/N0BIKS) 

BB2=20.*L0G<RHS2/RHS1)/106(10.) 

TYPE  10050,  R«S1,RHS2,BI2 
TYPE  10060 
ACCEPT  10070, GNDB 
TYPE  10075 
ACCEPT  10070, BINS 
TYPE  10080 
ACCEPT  10020, F3NM 
F3NAN<30)=0 

IF<.N0T.(F3NAN(1).NE,'  '))60  TO  9000 - >4BORT 


B-34 


PAGE  1 


C - ¥T*«<6  IT  MU  IICK 

C - friOH  VT)«6  29  ^JB  80 

C - MES  UEIGHTEB  MERGE  OF  TWO  FILES  AM  AMS  GAUSSIAN  WISE. 

SHE  F1MAN(30).F2NAN(30).F3NAM(30) 

B1NENSJQN  I5UF(254).JMJF(750)»KBUF<256) 

INTEGERS  Ik£C.I0SB»IMF<?)»JSBF<9bKBIF;9> 

EQUIVALENCE  <»gFiIMF).(JHF.JWF).(KMFiMUF) 

DATA  FlNAM/3OtO/F2NAN/3OIO/F3NAM/3Ot0/ 

DATA  KBUF/25410/ 

->  TfPE  10010 

ACCEPT  I0020.F1NAN 
F1NAN(30)=0 

OOPfN(UNIT=l*NAME=FlNAN»TYPEs'OU)'» 

1  8UFF£RC0UNT=-liERR=9000»REAB0NLY) 

ACCEPT  10020iF2NAN 
F2NAN(30)=0 

OOPEN<UNITs2»NANE=F2NAN»T'T£s'OLII,» 

1  BUFFERC0UNT=-1.ERR=9000.R£AMNLY> 

IREC=1 

CALL  REABF<1»IR£C,IBUF.512.I0SB) 

CALL  READF<2»IR£C.J8UF.512»I0SB) 

N1PLKS=IDBF(9) 

N2RLKS=JDBF<9) 

NSMKS-N1RLKS 

I F ( NSBLKS . 6T . N2BLKS)NSBLKS=N2BLKS 

TYPE  t0030.NSM.KS 

ACCEPT  10040.N0BLKS 

IF  ( NOELKS .  LE .  0  )NOMJ(S=NSBLKS 

IF<NOBL/<S.GT.NSRLKS>NOBLKS*NSBLKS 

RNS1=0. 

RNS2-0 . 

W  20  IBLK=1.N0BLKS 
h  IREC=I8LKT1 

CALL  READF(1,IREC,IBUF»512»I0SB) 

CALL  R£ADF<2.ItfC.JBUFi512rI0SB> 

BSU«1=0. 

BSUN2=0. 

DO  10  IDX-1.25A 
X^IBUFUDX) 

Bsum=BSummx 

X= JBOF ( IPX ) 

BSUN2=BSUH2+XtX 
10  CONTINUE 

RHSJsRHSHBSUHl/256. 

R«S2=RMS2fBSUM2/250, 

20  CONTINUE 

RNSl-SCRT(RftSl/NOBLKS) 

RHS2=SQRT ( RNS2/N0BLKS ) 

BB2r20 . *LOG<  RNS2/RNS1 )/LOG( 10 . ) 

TYPE  10050.  RNSl.WtS2.DB2 
TYPE  10060 
ACCEPT  10070.GK5B 
TYPE  '0075 
ACCEPT  10070.DBNS 
TYPf  10080 
ACCEPT  10020.F3MAM 
F3NAN<30)*0 

IF( .N0T.(F3NAN(1).NE.'  '))60  TO  9000 - >^BORT 
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WD«UMITi8.>W1E=F3W«»Trf>£='«tr » 

1  8UFFERCOUWT=-li£RR=9000) 

KMJFdM 

K3UF(2)=6400 

SEC=N0BLKS/25. 

MIN=N0B1KS/1500 

IS£C=SEC-60.*MIN 

K8UF<tt>=KIM 

KF4.IF(9)-ISFC 

KHhF<?)=fH3BLKS 

IREC=I 

call  NRITEF(8fIREC»KBUF*512»I0SB) 
RMSN=RMS1*10.»*!DBNS/20.) 

GAIN=10.tt!GNDB/20. ) 

FCTR1=1./!2.FGAIN> 

FCTR2=GAIN*FCTR1 
ISEED1M23 
ISEED?=454 
10  40  IBLK=1.N0BLKS 
A  IR£C=IBLK+1 

CALL  READF!l»IREC>IBUFi512»IOSB) 

CALL  READF(2>IREC,J3UF.512.I0S8) 

BO  30  IDX=1»256 
A  GSNS=0. 

DO  25  KGS=1»12 

f  GSNS=6SNSFRAN! ISEEDi » ISEED2) 

25  CONTINUE 

GSN5=RMSN*!CSNS-6.) 

IF!GSNS.GT. 32000. )GSN$=32000. 

IF (GSNS.LT. -32000. )6SWS=-32000. 

*  k  KBUF !  IDX  > =FCTR1 * ( IBUF ( IDX ) +GSNS ) +FCTR2I JBUF! IDX ) 

30  CONTINUE 

CALL  WRITEF!8,IREC,KBUF,512iI0SB) 

40  CONTINUE 

TYPE  10090/NQBLKS 

9000  CONTINUED - /ABORT 

•4-STOF 

100J0  FORMAT! '  TWO  FILES  FOR  WEIGHTED  MERGE  WITH  NOISE?') 
10020  FORMATdOAl) 

10030  FORMAT! '  HAVE'»I8»'  BLOCKS.  HCW  MANY  TO  MERGE?') 

10040  F0RMAKI8) 

10050  FORMAT!'  RMS'r2F10.2»'.  2WD/1STJ ' »F10.2» *  DB.) 

10060  FORMAT!'  WHAT  DB  GAIN  FOR  2ND  FILE?') 

10070  F0RMAT(F10.2) 

10075  FORMAT!'  WHAT  DB  6AUSS  NOISE  VS  1ST  FILE?(USE  .)') 
10080  FORMAT!'  WHAT  OUTPUT  FILE?') 

10090  FORMAT!'  HEADER  AND'fI8»'  DATA  BLOCKS  .‘RITTEN.') 

END 
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